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Abstract 

Divergence-conforming  B-splines  are  developed  for  application  to  the  incompressible 
Navier-Stokes  equations  on  geometrically  mapped  domains.  These  enable  smooth, 
pointwise  divergence- free  solutions  and  thus  satisfy  mass  conservation  in  the  strongest 
possible  sense.  Semi-discrete  methods  based  on  divergence-conforming  B-splines  are 
shown  to  conserve  linear  and  angular  momentum  and  satisfy  balance  laws  for  energy, 
vorticity,  enstrophy,  and  helicity.  These  are  geometric  structure-preserving  quantities 
and  numerical  simulations  that  are  sensitive  to  them  are  shown  to  be  qualitatively 
correct  and  quantitatively  accurate.  The  methods  developed  are  anticipated  to  open 
new  doors  to  the  practical  calculation  of  complex  flows  and  to  studies  of  their  physical 
behavior. 

Keywords:  Incompressible  Navier-Stokes  equations,  B-splines,  Isogeometric  analysis, 
Divergence-conforming  discretizations,  Structure-preserving  discretizations 


1  Introduction 

The  unsteady  incompressible  Navier-Stokes  equations  are  infused  with  vast  geometric 
structure,  evidenced  by  a  wide  array  of  balance  laws  for  momentum,  angular  momen¬ 
tum,  energy,  vorticity,  enstrophy,  and  helicity.  These  balance  laws  are  considered 
to  be  of  prime  importance  in  the  evolution  of  laminar  and  turbulent  flow  struc¬ 
tures  [38,  39,  40,  41],  and  they  are  even  believed  to  play  a  role  in  the  regularity  of 
Navier-Stokes  solutions  [8,  32],  The  key  to  unlocking  much  of  the  geometric  structure 
of  Navier-Stokes  flow  is  precisely  its  volume-preserving  nature,  yet  most  numerical 
methods  only  satisfy  the  incompressibility  constraint  in  an  approximate  sense.  Conse¬ 
quently,  such  methods  do  not  obey  many  fundamental  laws  of  physics.  In  particular, 
semi-discretizatations  which  conserve  momentum  are  typically  guaranteed  to  balance 
energy  if  and  only  if  the  incompressibility  constraint  is  satisfied  pointwise.  This  is 
especially  concerning  as  energy  plays  a  fundamental  role  in  numerical  stability  [35]. 
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In  this  paper,  we  present  new  divegence-conforming  B-spline  semi-discretizations 
for  the  unsteady  Navier-Stokes  problem.  These  semi-discretizations  are  motivated  by 
the  theory  of  isogeometric  discrete  differential  forms  [12,  13]  and  extend  the  steady 
Navier-Stokes  discretizations  presented  in  [23]  to  unsteady  Navier-Stokes  flows.  As 
incompressibility  is  satisfied  pointwise,  these  semi-discretizations  replicate  the  geo¬ 
metric  structure  of  the  unsteady  Navier-Stokes  equations  and  admit  discrete  balance 
laws  for  momentum,  angular  momentum,  energy,  vorticity,  enstrophy,  and  hclicity.  In 
this  sense,  our  semi-discretization  scheme  may  be  thought  of  as  a  structure-preserving 
or  mimetic  discretization  procedure  for  the  unsteady  Navier-Stokes  equations.  We  im¬ 
pose  no-penetration  boundary  conditions  strongly  and  no-slip  boundary  conditions 
weakly  using  Nitsche’s  method.  This  enables  our  method  to  handle  boundary  layers 
without  resorting  to  stretched  meshes  [6,  7].  This  also  allows  our  discretization  proce¬ 
dure  to  naturally  default  to  a  conforming  approximation  of  Euler  flow  in  the  limit  of 
vanishing  viscosity  and  to  possess  both  energy  and  helicity  as  inviscid  invariants.  The 
proposed  semi-discretizations  are  extended  to  general  mapped  geometries  of  scientific 
and  engineering  interest  using  divergence-  and  integral-preserving  transformations 
for  velocity  and  pressure  fields  respectively.  In  addition  to  all  the  features  mentioned 
above,  a  recent  paper  of  Guernrond  [24]  suggests  that  our  semi-discretizations  con¬ 
verge  to  physically  relevant  weak  solutions  satisfying  a  local  (in  space-time)  energy 
balance.  It  is  not  known  at  this  time  whether  or  not  such  a  convergence  property 
holds  for  spectral  semi-discretizations. 

The  use  of  B-splines  in  the  numerical  anlysis  of  unsteady  Navier-Stokes  flow  has 
already  been  conducted  with  much  success.  The  novelty  of  the  method  presented  here 
is  simply  the  use  of  tensor-product  B-splines  that  are  capable  of  exactly  satisfying  the 
incompressibility  constraint.  In  the  Direct  Numerical  Simulation  (DNS)  community, 
a  common  method  of  choice  in  simulating  wall-bounded  flows  is  the  use  of  Fourier 
spectral  discretizations  in  periodic  directions  and  B-splines  in  wall-normal  directions 
[33,  34,  37].  In  this  setting,  B-splines  are  often  preferred  over  polynomial-based  spec¬ 
tral  discretizations  due  to  their  high  resolving  power,  easy  implementation  of  bound¬ 
ary  conditions,  and  ability  to  employ  stretched  grids.  Recently,  Bazilevs  et  al.  studied 
the  turbulence  problem  in  a  series  of  papers  using  NURBS-based  isogeometric  analy¬ 
sis  in  conjunction  with  a  Variational  Multiscale  (VMS)  methodology.  In  these  papers, 
it  was  found  that  the  increased  continuity  of  splines  led  to  enhanced  numerical  results 
[1,  3,  4,  6,  7].  It  is  believed  that  much  of  this  success  can  be  attributed  to  the  spectral- 
like  properties  of  B-splines.  In  Figure  1,  we  have  plotted  the  phase  errors  associated 
with  one- dimensional  /c-method  NURBS  (which  in  this  setting  reduce  to  B-splincs 
of  maximal  continuity)  and  C°  finite  element  discretizations  of  the  first-order  wave 
equation.  Note  that  the  phase  error  associated  with  the  quadratic  NURBS  solution  is 
much  smaller  than  that  associated  with  the  quadratic  finite  element  solution.  Indeed, 
it  can  be  shown  that  the  phase  error  for  C P_1  NURBS  solutions  scales  like  0{h2p+ 2) 
while  the  phase  error  for  C°  finite  element  solutions  scales  like  0(h2p )  (see  Chapter  9 
of  [15]).  Recently,  the  theory  of  Kolmogorov  n- widths  has  been  utilized  to  shed  more 
light  on  the  approximation  properties  of  B-splincs  [20].  In  this  study,  it  was  revealed 
that  B-splines  are  much  more  accurate  on  a  per  degree-of- freedom  basis  than  stan¬ 
dard  finite  elements  and  possess  similar  approximation  properties  to  that  of  a  spectral 
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Figure  1:  The  first-order  wave  equation.  Phase  errors  versus  non-dimensional  wave 
numbers.  Comparison  of  linear  and  quadratic  finite  elements,  C 1  quadratic  NURBS, 
and  C2  cubic  NURBS. 


basis.  We  believe  that  by  combining  the  spectral-like  properties  of  B-splines  with  the 
preservation  of  the  geometric  structure  of  the  unsteady  Navier-Stokes  equations,  our 
semi-discretization  procedure  may  become  a  useful  tool  for  both  engineering  analysis 
and  the  mathematical  study  of  the  unsteady  Navier-Stokes  equations. 

An  outline  of  this  paper  is  as  follows.  In  the  following  section,  we  present  some 
basic  notation.  In  Section  3,  we  recall  the  unsteady  Navier-Stokes  problem  subject 
to  homogeneous  Dirichlet  boundary  conditions.  In  Section  4,  we  briefly  review  B- 
splines,  the  basic  building  blocks  of  our  new  discretization  technique,  and  in  Section 
5,  we  define  the  B-spline  spaces  which  we  will  utilize  to  discretize  velocity  and  pres¬ 
sure  fields.  In  Section  6,  we  present  our  semi-discrete  variational  formulation  for  the 
steady  Navier-Stokes  problem,  prove  well-posedness  and  a  discrete  energy  inequal¬ 
ity,  and  present  an  a  priori  error  estimate.  In  Section  7,  we  derive  various  balance 
laws  for  our  semi-discretizations,  including  balance  of  linear  and  angular  momentum, 
energy,  vorticity,  enstrophy,  and  helieity.  In  Section  8,  we  present  numerical  results 
illustrating  the  advantages  of  our  semi-discretization  procedure  on  three  benchmark 
problems:  two-dimensional  Taylor-Green  vortex  flow,  alternating  cylindrical  Couette 
flow,  and  three-dimensional  Taylor-Green  vortex  flow.  Each  of  these  problems  is  sen¬ 
sitive  to  preservation  of  conserved  quantities  and  the  growth  and  decay  of  functionals 
associated  with  geometrical  structure  of  the  flow.  In  Section  9,  we  draw  conclusions. 
Before  proceeding,  it  should  be  mentioned  that  we  do  not  consider  any  artificial 
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diffusion  mechanisms  or  subgrid  turbulence  models  in  this  paper.  As  such,  the  semi¬ 
discretizations  presented  here  should  only  be  utilized  if  the  flow  held  is  sufficiently 
resolved  by  the  spatial  mesh.  That  being  said,  standard  Large  Eddy  Simulation  mod¬ 
els  can  be  utilized  in  conjunction  with  the  proposed  semi-discretizations  to  capture 
fine-scale  turbulent  effects  on  coarse  meshes. 


2  Notation 


We  begin  this  paper  with  some  basic  notation.  For  d  a  positive  integer  representing 
dimension,  let  D  C  denote  an  arbitrary  bounded  Lipschitz  domain  with  boundary 
d D.  As  usual,  let  L2(D)  denote  the  space  of  square  integrable  functions  on  D  and  de¬ 
fine  L 2(D)  =  (L2(D))d.  We  denote  Lq(D)  C  L2(D)  as  the  space  of  square-integrable 
functions  with  zero  average  on  D.  We  will  also  utilize  the  more  general  Lebesgue 
spaces  LP(D )  where  1  <  p  <  oo.  Let  Hk(D)  denote  the  space  of  functions  in  L2(D) 
whose  kth- order  derivatives  belong  to  L2(D )  and  define  HA  (D)  =  (Hk(D))  .  We 
identify  with  Hk(D)  the  standard  Sobolev  norm 

(  ,  \ 1/2 

IMI/rfc(r>)  —  I  ^  ll-^avlli2(£>)  I 

\H<&  / 


where  ct  =  (an,  a2,  •  •  • ,  otd)  is  a  multi-index  of  non-negative  integers,  |a|  =  a±  +  ct2  + 
. . .  +  ad,  and 

Q\a\ 

dxl1  dx2 2  •  •  •  dxjd 

We  denote  the  Sobolev  semi-norms  as  |  •  | nk(D),  and  we  adopt  the  convention  H°(D)  = 
L2(D).  Throughout,  Sobolev  spaces  of  fractional  order  are  defined  using  function 
space  interpolation  (see,  e.g.,  Chapter  1  of  [48]).  We  define  H(\  (D)  c  Hl(D)  to  be 
the  subspace  of  functions  with  homogeneous  boundary  conditions  and  define  Hj(D) 
to  be  the  vectorial  counterpart  of  H(]  (D).  We  deffiie  H(div;  D)  to  be  the  Sobolev 
space  of  all  functions  in  ~L2(D)  whose  divergence  belongs  to  L2(D).  This  space  is 
equipped  with  the  norm 


|  H(div;D) 


=  V 


1 2 

L  2{D) 


+  l|diw|||2(jD) 


1/2 


We  also  define 


H0(div;  D)  =  {v  e  H(div;  D)  :  v  •  n  =  0  on  d D} 

where  n  denotes  the  outward  pointing  unit  normal.  Finally,  for  X  a  real  Banach 
space  and  S  a  positive  real  number,  we  define  Lq(0,  S;X)  as  the  space  consisting  of 
all  strongly  measurable  functions  (j)  :  (0,  S')  —±X  with 

/  S  1 

,s-,x)-=yJo  \\m\qxdt 
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for  1  <  q  <  oo  and 


||0IU-(O,S;X)  :=  ess  supo<t<s||0(t)||x  <  oo,  (2) 

and  C'°([0,  S];  X )  as  the  space  of  all  continuous  functions  (j>  :  [0,  S]  — >•  X. 


3  The  Unsteady  Navier-Stokes  Problem 


In  this  section,  we  recall  the  unsteady  Navier-Stokes  problem  subject  to  homogeneous 
Dirichlet  boundary  conditions.  For  d  a  positive  integer,  let  Q  denote  a  Lipschitz 
bounded  open  set  of  Throughout  this  paper,  d  will  be  either  2  or  3.  The  problem 
of  interest  is  as  follows. 


Given  v  G  R+,  f  :  SI  X  (0,  oo)  — »  Md,  and  u0  :  — >■  Md,  hnd  u  : 

0  x  [0,  oo)  — y  and  p  :  x  (0,  oo)  — »  M  such  that 


<9u  _  .  . 

—  +  V  •  (u  <g>  u) 


(S)l 


gradp  —  V  •  (2z/Vsu)  =  f 
divu  =  0 


in  £2  x  (0,  oo)  and 


u  =  0 

u(-,0)  =  u0(-) 


on  dQ  x  (0,  oo) 
in  Q. 


(3) 

(4) 

(5) 

(6) 


Above,  u  denotes  the  flow  velocity  of  a  fluid  moving  through  the  domain  Q,  p  denotes 
the  pressure  acting  on  the  fluid  divided  by  the  fluid  density,  v  denotes  the  kinematic 
viscosity  of  the  fluid,  f  denotes  a  body  force  acting  on  the  fluid  divided  by  the  density, 
and  Vsu  denotes  the  symmetrized  gradient  of  the  velocity  field  defined  by 

V*u  =  )  (Vu  +  (Vu)T)  . 


An  appropriate  definition  of  weak  solution  is  not  entirely  obvious  in  the  context 
of  unsteady  Navier-Stokes  flows,  especially  for  domains  in  M3.  The  most  basic  type 
of  weak  solution  is  a  so-called  Leray-Hopf  solution.  Given  a  fixed  end-time  T  >  0,  let 
us  assume  that  f  G  L°°(0,T;  LJ(h2))  and 


u0  G  {u  G  H0(div;  hi)  :  divu  =  0}  . 


A  Leray-Hopf  solution  over  the  time  interval  (0,  T)  is  then  defined  as  a  vector  function 
u  G  L°°(0,  T;  L2(fl))  fl  L2( 0,  T;  Hq(Q))  satisfying  u(-,  0)  =  u0,  divu  =  0  in  the  sense 
of  distributions,  the  energy  inequality 


u 


«(n)+  /  4H|Vsu(s)||^2(n))dxdGG  <  ||u(0)||22(n)+  /  2  (f(s), u(s))L2(n)  ds  (7) 
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for  almost  every  t  G  (0,  T),  and 


(u  •  dt-v  +  (u  (8)  u)  :  Vv  +  u  •  V  •  (2z/ VAv)  +  f  •  v)  (htdt  + 


/  u0  •  vdx 

Jfi 


0  (8) 


for  all  smooth  test  functions  v  G  (Cq°((0,  T)  xQ))d  such  that  divv  =  0.  In  1934,  Leray 
obtained  the  first  three-dimensional  existence  results  for  the  Navier-Stokes  equations 
by  carefully  constructing  a  Leray-Hopf  weak  solution  for  the  Navier-Stokes  problem 
posed  on  all  of  M3  [36],  and  seventeen  years  later,  Hopf  extended  this  existence  re¬ 
sult  to  the  Dirichlct  problem  by  constructing  a  Leray-Hopf  solution  as  the  limit  of 
a  sequence  of  Galerkin  approximations  [30].  While  the  question  of  existence  was 
answered  many  decades  ago  for  Leray-Hopf  solutions,  the  question  of  uniqueness  re¬ 
mains  unanswered  for  domains  in  M3.  This  is  primarily  due  to  an  intimate  relationship 
between  uniqueness  and  regularity  [18].  A  priori  estimates  have  not  yet  been  able  to 
preclude  the  occurrence  of  so-called  vorticity  bursts  reaching  scales  smaller  than  the 
Kolmogorov  scale  due  to  the  presence  of  enstrophy  production.  In  his  seminal  1977 
paper,  Scheffer  introduced  a  systematic  means  of  addressing  the  regularity  question 
by  studying  the  Hausdorff  measure  of  the  singular  set  of  weak  solutions  [45].  If  one 
can  prove  the  measure  of  this  set  is  zero,  one  will  then  have  answered  the  Navier- 
Stokes  smoothness  question  in  the  affirmative.  To  systematically  study  the  Hausdorff 
measure,  Scheffer  introduced  so-called  suitable  weak  solutions  which  satisfy  a  local 
space-time  energy  balance.  Specifically,  suitable  weak  solutions  satisfy  the  inequality 


dt  -|u|2  +  V-  Ulul2+P  u  -^V2  -|u|2  +HVu|2-f-u<0  (9) 


in  a  distributional  sense.  Such  a  local  balance  can  be  interpreted  as  an  “entropy 
condition”  for  incompressible  flows.  Suitable  weak  solutions  are  known  to  always 
exist,  and  they  may  be  obtained  via  a  regularization  of  the  Navier-Stokes  equations 
and  passing  to  the  limit.  With  the  notion  of  a  suitable  weak  solution,  Scheffer  was 
able  to  derive  a  bound  from  above  of  some  Hausdorff  measure  of  the  singular  set. 
This  bound  was  then  improved  upon  in  the  famous  1982  paper  of  Caffarelli,  Kohn, 
and  Nirenberg  where  it  is  proven  that  the  one-dimensional  Hausdorff  measure  of  the 
set  of  singularities  for  a  suitable  weak  solution  is  zero  [14].  That  is,  if  singularities 
exist,  they  cannot  lie  along  a  line  in  space-time.  This  is  widely  considered  to  be  the 
best  general  result  in  the  direction  of  the  Navier-Stokes  Millenium  Prize  Problem. 

All  of  the  above  machinery  is  not  necessary  in  the  derivation  of  a  semi-discrete 
variational  formulation.  As  we  will  show  later  in  this  paper,  we  are  able  to  obtain 
a  well-posed  semi-discrete  problem  via  a  simple  Galerkin  methodology.  That  being 
said,  a  recent  paper  of  Guernrond  suggests  that  our  semi-discrete  formulation  con¬ 
verges  to  a  suitable  weak  solution  [24],  It  should  be  noted  that  it  is  not  known 
whether  or  not  such  a  convergence  property  holds  for  spectral  discretizations  which 
are  the  standard  in  direct  numerical  simulation  of  turbulent  flows.  It  is  possible  that 
spectral  discretizations  may  not  converge  to  suitable  weak  solutions  as  they  do  not 
induce  enough  numerical  diffusion  to  counteract  the  Gibbs- Wilbraham  phenomenon, 
preventing  a  local  energy  balance  from  holding  in  the  limit  [25]. 
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4  B-splines  and  Geometrical  Mappings 

In  this  section,  we  briefly  introduce  B-splines,  the  primary  ingredient  in  our  discretiza¬ 
tion  technique  for  the  unsteady  Navier-Stokes  equations.  We  also  introduce  mappings 
which  will  allow  us  to  extend  our  discretization  technique  to  general  geometries  of 
engineering  interest.  For  an  overview  of  B-splines,  their  properties,  and  robust  algo¬ 
rithms  for  evaluating  their  values  and  derivatives,  see  de  Boor  [17]  and  Schumaker 
[46].  For  the  application  of  B-splines  to  finite  element  analysis,  see  Hollig  [29]  and 
Cottrell,  Hughes,  and  Bazilcvs  [15]. 

4.1  Univariate  B-splines 

For  two  positive  integers  k  and  n,  representing  degree  and  dimensionality  respectively, 
let  us  introduce  the  ordered  knot  vector 

S:={0  =  £i, (10) 


where 


Ci  6  £2  <  •  •  •  Cn+fc+1- 


Given  H  and  k,  univariate  B-splinc  basis  functions  are  constructed  recursively  starting 
with  piecewise  constants  (k  —  0): 


1  if  6  <  C  <  6+1 

0  otherwise. 


(ii) 


For  k  —  1,  2,  3, . . .,  they  are  defined  by  the  Cox-de  Boor  recursion  formula: 


£-6 
S i+k  Si 


+ 


Si+fc+l  C  r>k—l  / 

7  ~c  W+i 

Si+fc+l  Si+1 


(12) 


When  £i+k  —  =  0,  ^  is  taken  to  be  zero,  and  similarly,  when  6+fc+i  —  6+1  —  0, 

Ci+fc+i-g  js  taken  to  be  zero.  B-spline  basis  functions  are  piecewise  polynomials  of 

^ I  /c  |  1  h  1 

degree  k,  form  a  partition  of  unity,  have  local  support,  and  are  non-negative.  We 
refer  to  linear  combinations  of  B-spline  basis  functions  as  B-splines  or  simply  splines. 

Let  us  now  introduce  the  vector  C  =  {Ci;  •  •  •  >  Cm}  °f  knots  without  repetitions  and 
a  corresponding  vector  {rq, . . . ,  rm}  of  knot  multiplicities.  That  is,  rq  is  defined  to  be 
the  multiplicity  of  the  knot  Q  in  H.  We  assume  that  rt  <  k  +  1.  Let  us  further  assume 
throughout  that  r±  =  rm  —  k+1,  i.e.,  that  S  is  an  open  knot  vector.  This  allows 
us  to  easily  prescribe  Dirichlet  boundary  conditions.  At  the  point  Q,  B-spline  basis 
functions  have  oq  :=  k  —  r3  continuous  derivatives.  We  define  the  regularity  vector  ot 
by  a  :=  {0+, . . .  ,am}.  By  construction,  0+  =  am  =  —1.  In  what  follows,  we  utilize 
the  notation 

|d£  |  =  min{o;j  :  2  <  i  <  m  —  1}  (13) 

and  a  —  1  :=  (—1,  0.2  —  1, . . . ,  am_  1  —  1,  —1}  when  >  0  for  2  <  i  <  m  —  1. 
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We  denote  the  space  of  B-splines  spanned  by  the  basis  functions  £>f  as 

Sa  span  {Bj  }"=1  •  (14) 

When  k  >  1  and  ex*  >  0  for  2  <  %  <  m  —  1,  the  derivatives  of  functions  in  are 
splines  as  well.  In  fact,  we  have  the  stronger  relationship 

<i5) 

One  of  the  most  important  properties  of  univariate  B-splines  is  refinement  and,  per¬ 
haps  more  importantly,  nestedness  of  refinement.  Knot  insertion  and  degree  elevation 
algorithms  are  described  in  detail  in  Chapter  2  of  [15]. 


4.2  Multivariate  B-splines 

The  definition  of  multivariate  B-splines  follows  easily  through  a  tensor-product  con¬ 
struction.  For  d  again  a  positive  integer,  let  us  consider  the  unit  cube  O  =  (0,  l)d  C 
which  we  will  henceforth  refer  to  as  the  parametric  domain.  Mimicking  the  one¬ 
dimensional  case,  given  integers  ki  and  rq  for  /  =  1  let  us  introduce  open 

knot  vectors  H,  =  {£M, . . . ,  £,ni+kl+i,i}  and  the  associated  vectors  Ci  =  {(i,h  •  •  • ,  C 
{r1;/, . . . ,  and  ex;  =  {Cxi,;, . . . ,  There  is  a  parametric  Cartesian  mesh  M.h 

associated  with  these  knot  vectors  partitioning  the  parametric  domain  into  rectangu¬ 
lar  parallelepipeds.  Visually, 


M-h  —  {Q  —  <S>z=i,...,d  (( iuu  Cix+i.0  j  1  <  k  <  mi  ~  1}  • 


(16) 


For  each  element  Q  G  we  associate  a  parametric  mesh  size  h,Q  =  diarn(Q).  We 
also  define  a  shape  regularity  constant  A  which  satisfies  the  inequality 

A"1  <  <  A,  VQ  e  AV,  (17) 

h-Q 

where  hq  m jn  denotes  the  length  of  the  smallest  edge  of  Q-  A  sequence  of  parametric 
meshes  that  satisfy  the  above  inequality  for  an  identical  shape  regularity  constant  is 
said  to  be  quasi-uniform. 

We  associate  with  each  knot  vector  S;  (/  =  1 , ,d)  univariate  B-spline  basis 
functions  B^lt  of  degree  ki  for  %i  =  1, . . . ,  rq.  On  the  mesh  A4k,  we  define  the  tensor- 
product  B-spline  basis  functions  as 

Biu--£d  :=  Bh,i  ®  ®  Bid,d’  il  =  1’  •  •  •  ’  ni>  •  •  •  id  =  1’  •  •  •  ’  nd-  (18) 


We  then  accordingly  define  the  tensor-product  B-spline  space  as 


oki,...,kd 
<*  I,...,ad 


^\Z±d(Mh)  :=  span 


lki,...,kd 

,id 


ni,...,nd 

ii=l,...,id=l 


(19) 


Like  their  univariate  counterparts,  multivariate  B-spline  basis  functions  are  piece- 
wise  polynomial,  form  a  partition  of  unity,  have  local  support,  and  are  non-negative. 
Defining  the  regularity  constant 

a  :  =  min  min  {«*,;}  (20) 

1=1, ...,d  2<il<mi-l 

we  see  that  our  B-splines  are  C^-continuous  throughout  the  domain  Q.  Refinement  of 
multivariate  B-spline  bases  is  obtained  by  applying  knot  insertion  and  degree  elevation 
in  tensor-product  fashion.  In  the  remainder  of  this  paper,  we  consider  a  family  of 
nested  meshes  {-M-h}h<h0  an<^  associated  B-spline  spaces  {S^'"k^d(Mh)}h<ho  that 
have  been  obtained  by  successive  applications  of  knot  refinement.  Furthermore,  we 
assume  throughout  that  the  mesh  family  {-M-h}h<h0  Is  quasi-uniform. 

Note  that  each  element  Q  =  <S>z=i,...,d  (C ilti,  Ck+1,1)  has  the  equivalent  representation 
Q  =  ®i=i,...,d  (£ji,h  Cji+1,1)  f°r  some  index  jj.  With  this  in  mind,  we  associate  with  each 
element  a  support  extension  Q ,  defined  as 

Q  ■—  ®l=l,...,d{£ji—ki,h£ji+ki+l,l)  ■  (21) 

The  support  extension  is  the  interior  of  the  set  formed  by  the  union  of  the  supports 
of  all  B-spline  basis  functions  whose  support  intersects  Q.  Note  that  each  element 
belongs  to  the  support  extension  of  at  most  Hi=i,...,d(2ki  +  1)  elements. 


4.3  Piecewise  Smooth  Functions,  Geometrical  Mappings,  and 
Physical  Mesh  Entities 

On  the  parametric  mesh  A4h,  we  define  the  space  of  piecewise  smooth  functions  with 
interelement  regularity  given  by  the  vectors  alr..,  ol^  as 


_  /IOO 

,OLd  ^CXl,...,OLd 


(M, 


(22) 


Precisely,  a  function  in  ad  is  a  function  whose  restriction  to  an  element  Q  G  Mh 

admits  a  C°°  extension  in  the  closure  of  that  element  and  which  has  atl j  contin¬ 
uous  derivatives  with  respect  to  the  Zth  coordinate  along  the  internal  mesh  faces 
{(xi,...,xd)  :  xi  =  Ciuh  Cut'  <  xv  <  7^  0  for  all  k  =  2,..., mi  -  1  and 

ji 1  =  1, . . . ,  rri[i  —  1.  Note  immediately  that  any  function  lying  in  the  B-spline  space 

aiso  iies  in  cZ,..,ad- 

Unless  specified  otherwise,  we  assume  throughout  the  rest  of  the  paper  that  the 
physical  domain  O  C  l6*  can  be  exactly  parametrized  by  a  geometrical  mapping 
F  :  D  — »  Q  belonging  to  (C“  )  with  piecewise  smooth  inverse.  We  further 

assume  that  the  physical  domain  D  is  simply  connected  with  connected  boundary 
dfl  and  the  geometrical  mapping  is  independent  of  the  mesh  family  index  h.  The 
geometrical  mapping  F  naturally  induces  a  mesh 


JCh  =  {K:K=F(Q),QeMh} 


(23) 


on  the  physical  domain  D.  We  define  for  each  element  K  G  Kh  a  physical  mesh  size 


h'K  =  \\DF\\l°o(q)Jiq 


(24) 
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where  Q  is  the  pre-image  of  K,  and  we  also  define  the  support  extension  K  =  F(Q). 
We  define  for  a  given  mesh  the  global  mesh  size 

h  =  max  {hx,  K  e  /C^}  . 

Note  that  as  the  parametric  mesh  family  {M-h)h<ha  quasi-uniform  and  the  geomet¬ 
rical  mapping  F  is  independent  of  the  mesh  family  index  h,  the  physical  mesh  family 
{)Ch}h<ho  is  also  quasi- uniform.  We  refer  to  the  physical  domain  and  its  pre-image 

hi  interchangeably  as  the  patch.  It  should  be  noted  that,  in  general,  the  domain  hi 
cannot  be  represented  using  just  a  single  patch.  Instead,  multiple  patches  must  be 
employed.  For  the  sake  of  brevity,  the  multi-patch  setting  will  not  be  covered  in  this 
paper.  The  interested  reader  is  referred  to  [19,  22,  23] 

We  define  on  the  parametric  mesh  a  set  of  mesh  faces  Fh  =  {F}  where  F  is  a  face 
of  one  or  more  elements  in  JAh-  We  define  the  physical  set  of  mesh  faces  as 

Fh  =  {F=F(F):FeFh} 

and  we  define  the  boundary  mesh  to  be 

Th  =  {F  e  Th  :  F  C  912}  . 


By  construction, 


dtt  —  UperhF. 


Note  that  for  each  face  F  E  Th  there  is  a  unique  K  e  JCh  such  that  F  is  a  “face”  of 
K  (in  the  sense  that  F  is  the  image  of  a  face  of  Q,  the  pre-image  of  K).  We  hence 
define  for  such  a  face  the  mesh  size 


hp  :=  hx- 

One  may  also  define  hp  to  be  the  wall-normal  mesh-size  as  is  done  in  [6]. 

5  Discretization  of  Velocity  and  Pressure  Fields 

In  this  section,  we  define  the  B-spline  spaces  which  we  will  utilize  to  discretize  the 
velocity  and  pressure  fields  appearing  in  the  unsteady  Navier-Stokes  problem.  These 
spaces  are  motivated  by  the  recent  theory  of  isogeometric  discrete  differential  forms 
[12,  13]  and  may  be  interpreted  as  smooth  generalizations  of  Raviart-Thomas  elements 
[44],  For  a  more  in-depth  discussion  of  the  discrete  velocity  and  pressure  spaces  used 
in  this  paper,  see  Section  5  of  [22], 

5.1  Discrete  Spaces  on  the  Parametric  Domain 

Using  the  notation  of  the  previous  section  and  assuming  that 

a  :=  min{|cq|  :  l  =  1, . . . ,  d}  >  1, 
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we  define  the  following  two  spaces: 


V,  := 


,fc2 —  1 


fcl  — l,fc2 

a  2 


X  Q  1  ’ 

'ctl, a;2  —  l  ^  uai-l, 
qki,k2  —  1  ,/C3  —  1  w  Qfcl  —  I,fc2,fc3  —  1 

*-,cti,Q;2  — 1,CK3 


i  v  QKl  —  1,«2,B3  —  I  v  Qfcl-l,fe2-l,fc3 
-1  x  0cti-l,a2,«3-l  x  °ai-l,a2-l, 


OL  3 


if  d  —  2, 
if  d  —  3, 


Qh  := 


qkl  — 1,^2  —  1 

rf^i  — 1,A;2  — 1,^3  — 1 
^01—1,02  —  1,03  —  1 


if  d  —  2, 
if  d  =  3. 


The  spaces  V/,,  and  Qh  are  precisely  the  spaces  of  discrete  two-forms  and  three-forms 
that  were  introduced  in  [12],  and  they  were  first  applied  to  the  discretization  of 
Stokes  flow  in  [11].  The  space  V/,  comprises  our  set  of  discrete  velocity  fields  while  Qh 
comprises  our  set  of  discrete  pressure  fields.  Note  that  as  a  >  1,  our  discrete  velocity 
fields  are  H1 -conforming.  If  we  allow  a  =  0,  onr  spaces  collapse  to  standard  Raviart- 
Thomas  mixed  finite  elements  [44],  In  order  to  deal  with  no-penetration  boundary 
conditions,  we  make  use  of  the  following  constrained  discrete  spaces: 

Vo ,h  :=  { G  Vh  :  vh  ■  n  =  0  on  dQ 


Qo,h  ■=  1 5ft  e  Qh  :  J  qh(ht  =  0 

Above,  n  denotes  the  outward-facing  normal  to  dVt.  As  specified  in  the  introduction, 
we  choose  to  enforce  no-slip  boundary  conditions  weakly  using  Nitsche’s  method  [42], 
Due  to  the  special  relationship  given  by  (15),  the  spaces  Vo y  and  Q0,h  along  with  the 
parametric  divergence  operator  form  the  bounded  discrete  cochain  complex 

p;  div  pi 

Vo  ,h  - >  WiO.h 

where  div  is  the  divergence  operator  on  the  unit  cube  Q. 


5.2  Discrete  Spaces  on  the  Physical  Domain 

To  define  onr  discrete  velocity  and  pressure  spaces  on  the  physical  domain,  we  intro¬ 
duce  the  following  pullback  operators: 


t„(v)  :=  det(DF)  (DF)-1  (voF) ,  veH„(div;f!) 

ip[q)  ■-  det  (DF)  (qoF) ,  q  €  L%(S2) 


(25) 

(26) 


where  DF  is  the  Jacobian  matrix  of  the  parametric  mapping  F.  With  these  operators 
defined,  we  have  the  following  commuting  diagram: 


H0(div;  Q) 

Lit, 

H0(div;  Q) 


div 


♦  Lm 


(27) 


div 


*  Lim. 


ii 


This  motivates  the  use  of  the  following  discrete  velocity  and  pressure  spaces  in  the 
physical  domain: 

V0 ,h  ■=  jv  6  H0(div;  hi)  :  tu(v)  G  V0,/,j  , 

Qo,h  |<7  £  L20(n)  :  ip(q )  G  Qo,/,  j  • 

In  [12],  it  was  shown  that  there  exist  projection  operators  Ily  :  H0(div;  Q)  — >  Vo,/, 
and  IIgh  :  L2(Q)  — >  Q0,/i  such  that  the  following  proposition  holds. 

Proposition  5.1.  The  following  diagram  commutes: 


H0(div;  Fl) 


div 


»  Lim 


n° 


V, 


div 


0  ,h 


n°Q, 


>  Q 


(28) 


0,h- 


Furthermore,  there  exists  a  positive  constant  Cu  independent  of  h  such  that 
^vIIhPo)  <  C'nllvIlHho),  Vv  G  H0(div;  FI)  n 


in? 


(29) 


We  immediately  have  an  inf-sup  condition  for  our  discrete  velocity/pressure  pair 
as  a  consequence  of  the  above  proposition.  A  complete  proof  of  the  inf-sup  condition 
may  be  found  in  [22], 

Proposition  5.2.  There  exists  a  positive  constant  j3  independent  of  h  such  that: 

inf 

qh£Qo,h 

We  also  have  the  following  result. 

Proposition  5.3.  If  v/,  G  Vo,/,  satisfies 

(divvft,  qh)L2(p.)  =  0,  Wqh  G  Q0}h,  (31) 


sup 

vhSVo,ft 


(divv/,.,  qh)L2(n) 
lv^llH1(n)lkft.|U2(n) 


>P- 


(30) 


then  divv/,  =  0. 


Proof  The  proof  holds  trivially  as  div  maps  Vo./,  onto  Qo,h-  □ 

Hence,  by  choosing  Vo,/,  and  Qo,/,  as  discrete  velocity  and  pressure  spaces,  we  arrive 
at  an  inf-sup  stable  discretization  that  automatically  returns  velocity  fields  that  are 
pointwise  divergence- free. 
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5.3  Approximation  Results  and  Trace  Inequalities 

Let  us  define 

k' =  min  \ki  —  1\ .  (32) 

Note  that  the  discrete  velocity  and  pressure  spaces  Voy  and  Q0,/i  consist  of  mapped 
piecewise  polynomials  which  are  complete  up  to  degree  k! .  The  following  result  details 
the  local  approximation  properties  of  our  discrete  spaces.  Its  proof  may  be  found  in 
[12], 

Proposition  5.4.  Let  K  £  /Q,  and  K  denote  the  support  extension  of  K .  For 
0  <  j  <  s  <  k!  +  1,  we  have 

|v  -  n“,  v|HJ(K)  <  Ch^l|v||H.(*),  Vv  e  H'(A')  n  H„(div;  SI)  (33) 

li-na.«l»m<ChrikllH.(^,  Vq£H-(K)nLl(Sl)  (34) 

where  C  denotes  a  positive  constant,  possibly  different  at  each  appearance,  independent 

°f  h. 

Hence,  our  discrete  spaces  deliver  optimal  rates  of  convergence  from  an  approxi¬ 
mation  point  of  view.  We  will  also  need  the  following  trace  estimate  in  what  follows. 
Its  proof  can  be  found  in  [19]. 

Proposition  5.5.  Let  K  £  /Q.  Then  we  have 

II  (Vsvh)  n||  {L2{dK))d  <  Ctracehff  1 1  vh  1 1  h1  (K) ,  VvA  £  V0)h  (35) 

where  Ctrace  denotes  a  positive  constant  independent  of  h. 

In  [21],  it  was  shown  that  Proposition  5.5  holds  for  B-splines  and  parametric  finite 
elements  with  Ctrace  ~  (k')2.  However,  our  numerical  experience  has  indicated  that 
a  corresponding  global  trace  inequality  holds  with  Ctrace  ~  k!  if  B-splines  of  maximal 
continuity  are  utilized. 

6  The  Semi-Discrete  Problem 

In  this  section,  we  approximate  the  unsteady  Navier-Stokes  problem  using  the  discrete 
velocity  and  pressure  spaces  introduced  in  the  previous  section.  We  prove  the  resulting 
semi-discretization  scheme  is  well-posed  and  satisfies  a  discrete  energy  balance  law, 
and  we  briefly  discuss  a  priori  error  estimates. 

6.1  Semi-Discrete  Variational  Formulation 

We  begin  this  section  by  presenting  a  semi-discrete  variational  formulation  for  the 
unsteady  Navier-Stokes  equations  subject  to  homogeneous  Dirichlet  boundary  condi¬ 
tions.  Since  members  of  Vo,h  do  not  satisfy  homogeneous  tangential  Dirichlet  bound¬ 
ary  conditions,  we  employ  Nitsche’s  method  to  weakly  enforce  no-slip  boundary  con¬ 
ditions.  To  this  effect,  we  define  the  following  bilinear  form  for  Cpen  >  1  a  chosen 
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positive  penalty  constant  independent  of  mesh  size: 


h( w,  v)  =  k( w,  y)-Y,  /  2 v 

fg rh  Jf 

where 

k( w,v)  =  (2iA7sw,  Vsv)(Z;2(n))dxd ,  Vw,  v  G  H^Q). 
Let  us  also  define  our  discrete  space-time  velocity  space  as 


((Vsv)  n)  •  w+  ((Vsw)  n) 


a 


•  v  — 


pen 


hr 


w  •  V 


ds 

(36) 

(37) 


Vt  -  {v/t  €  C°((0,T);Vo,»)  :  d,vh  €  i2((0,  T);  V„,„)}  (38) 

and  our  discrete  space-time  pressure  space  as 


ay=  i2((o,r);so,„). 


(39) 


It  remains  to  specify  our  discrete  initial  condition.  In  what  follows,  we  choose 
Uo^  =  Ily  u0.  Note  that,  by  construction,  div u0:/,,  =  0.  We  can  alternatively  specify 
our  initial  condition  through  H1-projection  into  the  B-spline  space  of  divergence-free 
fields.  This  amounts  to  solving  a  steady  Stokes  problem  as  a  pre-processing  step. 
Employing  the  terminology  defined  above,  our  semi-discrete  formulation  reads  as  fol¬ 
lows. 

Find  u h  €  and  ph  G  Q ?  such  that  u^(0)  =  u0i/,  and,  for  almost 
every  t  G  (0,T), 


(G)< 


(dtnh(t),  vA)L2(n)  +  kh(uh(t),vh) 

+c(uh(t),uh(t)-,vh)  —  b(ph(t) ,  v h)  +  b(qh,uh(t))  =  (f(t),  vh)La(n) 

(40) 

for  all  Vh  £  Vo,/i  and  qh  G  Qo,h  where 

b(q,v)  =  (q,  divv)L2(Q) ,  Vg  G  L2(fl),  v  G  H^fl) ,  (41) 


c(w,  x;  v)  =  —  (w  ®  x,  Vv) 


(L2(U))dxd 


,  Vw,  x,  v  G  H1(fi) . 


(42) 


We  immediately  remark  that  the  semi-discrete  problem  we  have  obtained  is  a  set 
of  coupled  first-order  nonlinear  ordinary  differential  equations.  As  such,  we  can  use 
standard  approaches  from  the  theory  of  ordinary  differential  equations  to  obtain  exis¬ 
tence,  uniqueness,  and  regularity  results.  We  may  also  use  standard  time-integration 
schemes  to  obtain  a  fully-discrete  formulation. 

We  have  the  following  lemma  detailing  the  consistency  of  our  numerical  method 
for  sufficiently  regular  exact  solutions. 


Lemma  6.1.  Suppose  that  a  Leray- H op f  solution  (u ,p)  of  the  homogeneous  unsteady 
Navier-Stokes  problem  satisfies  the  regularity  conditions 

dtu  G  L2(0,  T;  L2(fl)),  u  G  L2(0,  T;  H3/2+e(fi)) 
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where  e  >  0.  Then,  for  almost  every  t  G  (0,  T): 

(0tu(t),  vh)L2(n)  +  kh(u(t),  vA)  +  c(u(t),  u(f);  vft) 

-6(p(<),vh)  +  6(gA,u(i))  =  (f(i),vh)L2(n} 

(43) 

for  all  \h  G  V0,h  and  qh  G  Q0,h- 

Proof.  We  trivially  have,  for  almost  every  t  G  (0,  T), 

&(gh,u(t))  =  0,  G  Q0)h.  (44) 

Now  let  Wh  G  Vo, h-  By  the  Sobolev  trace  theorem,  our  smoothness  assumption  for  u 
guarantees  that  (Vsu)  n  is  well-defined  along  <90  and  belongs  to  L2(0,  T;  (L2(<9fi))d). 
Hence,  we  can  utilize  integration  by  parts  and  the  fact  that  u(f)  satisfies  homoge¬ 
neous  Dirichlet  boundary  conditions  and  satisfies  homogeneous  normal  Dirichlet 
boundary  conditions  to  write,  for  almost  every  t  G  (0,T), 

(dtu(t),  vh)L2(Q)  +  kh(u(t),vh)  +  c(u(t),  u (t);  vh)  -  b(p(t),vh) 

=  f  (dtu(t)  -  V  •  (2z/Vsu(t))  +  V  •  (u  (8)  u)  +  Vp(t))  ■  wh 
Jn 

=  (f  W,v42(fi)  • 

This  completes  the  proof  of  the  lemma.  □ 

We  also  have  the  following  lemma  which  is  a  direct  consequence  of  Proposition 
5.3. 

Lemma  6.2  (Conservation  of  mass).  Suppose  (u h,Ph)  G  Vf  x  Qf  is  a  solution  to 
( G ).  Then 

\\divuh(t)\\2L2(n)dt  =  0.  (45) 

That  is,  divu/(  =  0  as  a  distribution. 

Weak  imposition  of  no-slip  boundary  conditions  allows  our  methodology  to  au¬ 
tomatically  default  to  a  compatible  semi-discretization  of  Euler  flow  in  the  setting 
of  vanishing  viscosity.  Moreover,  for  large  Reynolds  number  flows,  there  is  a  sharp 
boundary  layer  in  the  vicinity  of  walls.  Utilzing  Nitsche’s  method  allows  us  to  ac¬ 
count  for  these  layers  in  a  stable  and  consistent  manner  without  having  to  directly 
resolve  them  [5,  6,  7].  In  fact,  Nitsche’s  method  can  be  interpreted  as  a  variationally 
consistent  wall  model.  To  better  see  this,  let  us  formally  rewrite  our  semi-discrete 
variational  equations  as 

(dtuh(t),  vfc)L2(n)  +  f  T(t)  :  Vavhdx-  V  [  Q(t)-vhds 
Jn  FeFh  Jf 

+c(uh(t),  uft(t);  vft)  -  b(ph(t),vh )  +  b(qh,  u h(t))  =  (f(t),  va)l2(0)  (46) 
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where  T  is  a  symmetric  tensor  satisfying 

[  T(t)  :  Wdx  =  /  2z/Vsu h{t)  :  Wdx  - 


=  /  2z/u/j(t)  •  divWdx 


^  /  2z/u/l(t)  •  (Wn)  ds 
Ferh  -'F 

(in  the  sense  of  distributions) 


(47) 


Jn 

for  symmetric  tensors  W  with  well-defined  normal  trace  and  Q  is  a  vector  satisfying 

Q (t)  =  ^(Vsu fit))  n  -  ^^ufit)^  .  (48) 

Above,  T  is  a  weakly  defined  viscous  stress  tensor  and  Q  is  the  resultant  viscous 
boundary  traction  vector.  The  tangential  component  of  Q,  given  by 


Q tang  =  Q  -  (Q  •  n)  n, 


(49) 


is  the  effective  wall  shear  stress  vector.  As  the  semi-discrete  velocity  field  satisfies  the 
no-penetration  boundary  condition  strongly,  the  vector  Q  is  equal  to  the  semi-discrete 
shear  stress  2v  (Vsu^)  n  plus  an  additional  wall  shear  stress  term  Q+  in  the  direction 
tangent  to  the  wall.  Specifically,  we  have 


where 


=  — u 


*2 


u*2  = 


ypen 


h 


U  h 

1  II 

(50) 

Ml 

IIMI 

(51) 

For  under-resolved  flow  simulations,  the  magnitude  of  (V*^)  n  in  the  direction  tan¬ 
gent  to  the  wall  is  relatively  small  and,  as  such,  the  tangential  component  of  Q  is 
dominated  by  Q+.  In  this  sense,  Q+  becomes  a  model  for  the  wall  shear  stress.  As 
the  mesh  is  refined  and  the  flow  is  resolved,  Q+  — >  0.  The  above  interpretation 
allows  us  to  design  physically  motivated  penalty  values  for  Nitsche’s  penalty  param¬ 
eter.  Notably,  u*  may  be  interpreted  as  the  friction  velocity.  By  specifying  the  value 
of  u*  using  Spalding’s  law  of  the  wall  [47],  we  recover  a  standard  wall  model  for 
under- resolved  flow  simulations.  For  more  on  this  approach,  see  Section  3  of  [6].  It 
is  important  to  note,  however,  that  the  numerically  inspired  version,  (50),  produced 
results  of  the  same  quality  as  the  u*  given  by  Spalding’s  physically  inspired  law  of 
the  wall. 


Remark  6.1.  If  we  wish  to  impose  non-homogeneous  tangential  Dirichlet  ( e.g pre¬ 
scribed  slip)  boundary  conditions,  we  must  add  the  following  expression  to  the  right 
hand  side  of  our  semi- discrete  formulation: 

-  £  /F2" 

fg rh  JF 

where  u Bc  is  a  prescribed  vector  function  defined  on  <912  with  prescribed  boundary 
values.  If  we  also  wish  to  impose  non-homogeneous  normal  Dirichlet  (e.g.,  prescribed 


-  ((Vsvfe)  n)  •  uBC(t)  +  Bc(t)  ■  vh  )  ds  (52) 

it  pr 
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penetration)  boundary  conditions,  we  must  impose  these  strongly  and  add  the  following 
expression  to  the  left  hand  side  of  our  semi-discrete  formulation: 


CuwiUH{t),Vh)  =  Y,  /  (ubcO)  •  n 

F£ rh  Jf 


u h{t)  ■  vhds 


(53) 


and  the  following  expression  to  the  right  hand  side  of  our  semi-discrete  formulation: 


where 


and 


fuw(yh )  =  -  ^2  /  (uBc(t)  ■  n)_  uBC(t)  ■  vhds 
Fe rh  Jf 


(u BC(t)  ■  n) ,  = 


(u Bc(t)  ■  n)_  = 


ubcW  •  n  if  uBc(t)  •  n  >  0 


0 


otherwise 


u Bc(t)  ■  n  if  uBc(t)  •  n  <  0 


0 


otherwise. 


These  additional  terms  correspond  to  upwinding. 


(54) 


6.2  Energy  Balance 


Our  semi-discrete  formulation  satisfies  a  discrete  energy  balance  law  if  we  assume  that 
the  constant  in  Nitsche’s  method  is  chosen  large  enough.  First,  in  light  of  Proposition 
5.5,  we  assume  that 


(-  'pen  —  PoinC Korn ' 


(V6Vfe)nj \\L2{dK)y 


vJ 


VK  G  JCh,  vh  G  V, 


0  ,h 


(55) 


where  CpQin  is  the  positive  constant  appearing  in  Poincare’s  inequality: 


IMInqn)  <  CPoin\v |Hi(n),  Vv  G  H1(Q)  n  H0(div;O) 

and  C Korn  is  the  positive  constant  associated  with  the  following  Korn’s  inequality 

[10]: 

I  Wl  H1  (0)  <  CKorn  ( ||' W  ||  (L2(n))dxd  +  |' ||  w||JL2(an))d)  ,  Vw  G  H^ft). 


Second,  we  assume  that 

Cpen>Ahv\d(l\-lW-1')  (56) 

where  ho  is  the  mesh  size  of  the  coarsest  mesh  /Co  and  <90  denotes  the  length  of  d(l 
for  d  =  2  and  the  area  of  d(l  for  d  —  3.  Assumptions  (55)  and  (56)  ensure  that  the 
bilinear  form  kh(-,  •)  is  coercive  (see  [22,  23]  for  details). 

Lemma  6.3  (Global  balance  law  for  energy).  Suppose  (u  h,Ph)  £  VyX  Qff  is  a  solution 
to  ( G ).  Furthermore,  assume  (55)  and  (56)  are  satisfied.  Then 


1  d 

2  dt 


uh(f)llL2(Q) 


-kh(uh(t),  u h(t))  +  (f(t),  ufc(i))L2(n)  <  (f(t),  uh(i))L2(n)  (57) 


for  almost  every  t  G  (0,  T). 
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Proof.  Insert  (vh,Qh)  =  (uh(f),0)  into  (40)  for  almost  every  t  G  (0,  T)  and  use  the 
fact  that  clivus  =  0  as  a  distribution  to  obtain  the  expression 

(<9tu  h(t),uh(t))L2(n)  +kh(uh(t),uh(t))  +  c(uh(t),uh(t);uh(t))  =  (f(t),uh(t))L2(a). 

An  immediate  application  of  the  product  rule  yields 

^II^WIIl2^)  +  kh(uh(t),uh(t))  +  c(uh(t),uh(t);  uh(t))  =  (f(t),uh(t))L2(Q)  . 

To  proceed,  we  write 


c(uh(t),  u h{t)]  u h{t))  =  -  (u h{t)  <g>  u h{t))  :  Vu h(t)dx. 

Jn 

Since  divu^  =  0,  we  have 

c(uh(t),uh(t);uA(t))  =  I  div  (uh(t)|uft(t)|2)  dx, 

*  JQ 

and  by  the  divergence  theorem, 


c(uh(t),uh(t);uh(t))  =  0. 

To  complete  the  proof,  note  that  kh(-,  •)  is  coercive  due  to  Corollary  6.2  of  [23].  □ 

The  above  energy  balance  law  is  analogous  to  the  balance  law  satisfied  by  Leray- 
Hopf  weak  solutions.  However,  while  we  have  an  equality  in  our  discrete  balance 
law,  the  balance  law  for  Leray-Hopf  weak  solutions  is  characterized  by  an  inequality. 
This  inequality  is  an  artifact  of  regularization  and  passing  to  the  limit.  Note  that 
when  the  applied  forcing  term  is  conservative  (i.e.,  f  =  Vq  for  some  scalar  potential 
q  :  — >  R),  we  have 

^K«llh(n)  <  o 

for  almost  every  t  G  (0,  T).  Hence,  our  formulation  properly  dissipates  energy  in  the 
presence  of  conservative  forces.  Furthermore,  when  the  viscosity  is  taken  to  vanish, 
we  obtain 

Thus,  just  as  in  the  infinite- dimensional  setting,  energy  is  an  inviscid  invariant  for  our 
semi-discrete  formulation.  Consequently,  our  formulation  exhibits  time-reversibility. 


6.3  Global  Existence  and  Uniqueness 

Using  standard  approaches  from  the  theory  of  ordinary  differential  equations,  we 
obtain  the  following  theorem. 

Theorem  6.1.  Assume  (55)  and  (56)  are  satisfied.  Then  Problem  ( G )  has  a  unique 
solution  (u hiVh)  G  V!f  x  Qfi.  Moreover, 

II U^i  II  L°°(0,T;L2(r2))  ^  ^  ( II U0,h  ||  L2(f2)  +  l|f  II L2(0,T;L2(O)))  ■  (58) 
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Proof.  We  begin  by  establishing  existence  and  uniqueness  for  the  semi-discrete  veloc¬ 
ity  solution.  To  do  so,  we  restrict  ourselves  to  the  divergence-free  space 

Vj’  :=  {vh  G  :  divv/i  =  0}  . 

Each  function  G  Vf  can  be  written  uniquely  as 

m 

vh(t)  =  y^ej(t)wj 

i=  1 

where  {w,}^  is  a  basis  for  the  space 

V0 ifl  :=  {w h  G  V0,h  ■  divwA  =  0}  . 

Representing  our  desired  solution  as 


m 

nh(t)  =  Wj, 

i= 1 

the  semi-discrete  problem  in  the  kernel  becomes  a  nonlinear  system  of  first-order  or¬ 
dinary  differential  equations  for  the  coefficients  d*(f)  subject  to  appropriately  defined 
initial  conditions.  The  energy  estimate  in  Lemma  6.3  can  then  be  used  in  conjunction 
with  Gronwall’s  inequality  to  show  that  there  exists  a  unique  absolutely  continuous 
function  d(t)  =  (di(t), ,  dm(t ))  such  that  the  initial  conditions  are  satisfied  and  the 
nonlinear  system  of  first-order  ordinary  differential  equations  is  satisfied  for  almost 
every  time  t  G  (0,  T).  Hence,  a  solution  G  Vf  exists  and  is  unique.  Existence  and 
uniqueness  of  p\x  G  Off  is  an  immediate  consequence  of  the  inf-sup  condition  for  the 
bilinear  form  &(•,  •).  To  obtain  the  L 2  stability  estimate,  we  write  using  Lemma  6.3 

^KWIIlRq)  <2(f(t),uh(i))La(n)  <  ||f(*)||L2(n)  +  IM^IIl^o) 

for  almost  every  t  G  (0,  T).  The  desired  estimate  is  then  an  immediate  result  of  the 
differential  form  of  Gronwall’s  inequality.  □ 

Note  that  the  above  theorem  gives  existence  and  uniqueness  results  for  any  end- 
time  T.  Hence,  we  have  global-in-time  existence  and  uniqueness  for  our  semi-discrete 
problem.  Moreover,  our  method  is  well-posed  in  the  sense  that  it  returns  semi-discrete 
solutions  which  depend  continuously  on  the  given  data.  This  gives  our  methodology 
firm  mathematical  grounding. 

6.4  A  Priori  Error  Estimates 

If  one  assumes  that  a  Leray-Hopf  solution  to  the  unsteady  Navier-Stokes  equations 
is  sufficiently  smooth,  one  can  use  standard  functional  analysis  techniques  to  prove 
our  semi-discrete  velocity  solution  converges  at  optimal  order  (with  respect  to  the 
mesh  size)  to  the  Leray-Hopf  velocity  solution.  The  main  idea  is  due  to  Heywood  and 
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Rannacher  [26]  and  involves  splitting  the  velocity  error  into  two  components:  an  error 
contribution  associated  with  the  semi-discrete  solution  of  a  linearized  unsteady  Stokes 
problem  and  a  remainder  term  comprising  the  difference  between  the  semi-discrete 
Stokes  solution  and  the  semi-discrete  Navier-Stokes  solution.  The  resulting  analysis 
is  lengthy,  tedious,  and  tangential,  and,  as  such,  we  have  elected  not  to  include  it 
here.  The  final  a  priori  estimate  is  reported  in  the  theorem  below,  and  the  interested 
reader  is  directed  to  Chapter  9  of  [19]. 

Theorem  6.2.  Assume  that  (55)  and  (56)  are  satisfied  and  that  the  domain  hi  satisfies 
the  elliptic  regularity  condition.  Let  u  denote  a  Leray-Hopf  solution  of  the  unsteady 
Navier-Stokes  equations  and  ( uh,ph )  denote  the  unique  weak  solution  of  ( G ).  If  the 
regularity  conditions 

u  G  L°°(0,  T ;  HJ+1(0)),  dt u  G  L2(0,T;  HJ+1(0)) 
hold  for  j  >  1,  then  we  have 

llU  UCIl°°(0,T;L2(Q))  — 

cNS(i  +  <3Texp{/?r})/i2(*+1>  (|MI£„(0iT.h.+i(!!))  +  (2^)-1||a,u||2a(0T.H,+.(n))) 

(59) 

where  s  =  min  {k',j},  Q  and  [3  are  positive  constants  independent  of  h  and  T  that 
scale  asymptotically  with  z/-1  and  depend  on  the  supremum  of  u  over  kl  x  [0,  T),  and 
Cns  is  a  positive  constant  independent  of  h,  v,  and  T  that  scales  asymptotically  with 
r 

w pen  • 

Note  that  while  the  above  theorem  gives  an  optimal  L2  error  bound  for  the  velocity 
held  in  terms  of  the  mesh  size  h,  the  bound  grows  exponentially  in  time.  In  general, 
such  a  dependence  is  unavoidable.  However,  if  we  assume  an  exact  solution  is  stable 
in  the  sense  that  perturbations  decay  exponentially  in  time,  we  can  use  standard 
techniques  to  prove  error  estimates  which  are  uniform  in  time  [27].  Furthermore,  it  is 
known  that  solutions  of  the  unsteady  Navier-Stokes  problem  may  experience  blow  up 
at  the  initial  time  unless  certain  compatibility  conditions  relating  the  initial  condition 
and  applied  forcing  are  satisfied.  The  inherent  smoothing  properties  of  unsteady 
Navier-Stokes  flow  may  be  exploited  in  order  to  derive  optimal  error  estimates  away 
from  the  initial  time  for  higher-order  spatial  discretizations  [28].  Finally,  an  additional 
analysis  can  be  conducted  to  obtain  error  estimates  which  are  suboptimal,  by  one 
order,  for  the  semi-discrete  pressure  held.  Our  numerical  experience  has  suggested 
that  the  semi-discrete  pressure  held  converges  at  optimal  order  in  contrast  with  these 
estimates. 


7  Balance  Laws 

In  this  section,  we  present  a  collection  of  balance  laws  for  our  semi-discretization 
scheme.  These  balance  laws  supplement  the  discrete  energy  balance  law  derived  in 
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Subsection  6.2  and  give  our  semi-discrete  formulation  geometric  structure.  Before 
proceeding  further,  we  pose  our  semi-discrete  problem  in  a  slightly  different  form. 
This  form  will  more  directly  reveal  the  conservation  structure  of  our  chosen  semi¬ 
discretization.  Let  us  introduce  the  following  discrete  trace  space: 

Vtrace,h  ■=  {<?  e  L2  (dQ)  :  q  =  vh  ■  n,  vh  G  V/(}  .  (60) 


Note  that,  in  the  above  expression,  Vh  denotes  the  space  of  discrete  velocity  fields 
which  are  free  of  Dirichlct  boundary  conditions.  Vtrace,h  represents  the  natural  trace 
space  associated  with  Vh,  and  it  is  a  Hilbert  space  when  endowed  with  the  standard 
L 2  inner-product  over  dil.  Introducing  the  space-time  trace  space 

Vtace,T  :=  L2(0,T-Vtrace,h),  (61) 


let  us  consider  the  following  semi-discrete  problem. 

'  Find  uh  G  V£,  ph  G  Q £,  and  ptrace,h  e  VthraceT  such  that  uh(0)  =  u0)h 
and,  for  almost  every  t  G  (0,  T), 

(dtuh(t),vh)L2{n)  +  kh(uh(t),vh)  +  c(uft(t),uft(t);  vft) 

( v4.')  < 

-b(ph(t),Vh)  +  (ptraee,h(t),Vh-n)L2{dQ)  +  b{qh,Uh(t))  =  (f(t),  Vfc)La(n) 

(62) 

,  for  all  \h  e  Vh  and  qh  G  Q0,h- 

Here,  we  have  introduced  the  auxiliary  field  ptrace,h  and  released  the  no-penetration 
boundary  condition  on  the  discrete  test  space  for  the  momentum  equations.  Note  that 
the  solution  of  Problem  (G)  is  also  a  solution  to  Problem  (H)  (modulo  the  auxiliary 
held  Ptrace,h )  since 


i,Ptrace,h{t')  ,^h  '  ^)_t2(gn)  —  0)  G  Vo ,h- 

Furthermore,  the  auxiliary  field  ptrace,h  is  unique  due  to  the  obvious  inf-sup  condition 

inf  sup  {qh,wh  ■  n)L2(9Q)  =  |M|£2(sn)- 

trace, h 

By  employing  integration  by  parts,  we  observe  the  auxiliary  field  ptraCe,h  approximates 
the  trace  of  the  pressure  field.  Hence,  Ptrace,hn  gives  the  semi-discrete  normal  traction 
on  <9h2  due  to  pressure  forces. 


7.1  Conservation  Properties  on  Rectilinear  Domains 

Suppose  that  H  is  a  rectilinear  domain  that  has  been  mapped  from  parametric  space 
using  an  affine  transformation.  Then,  the  unit  vectors  e,:  G  necessarily  belong 
to  the  discrete  space  Vh-  If  we  select  wh  =  e*  and  qh  =  Phif)  in  (62)  and  sum  over 
%  —  1, . . . ,  d,  we  obtain  the  following  discrete  balance  law  for  linear  momentum. 
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Global  conservation  of  linear  momentum  on  rectilinear  domains.  For  recti¬ 
linear  domains,  the  semi- discrete  velocity  field  satisfies 

4  f  u h(t)dx.=  I  (-Ptrace,h(t)n  +  Q(t))  ds  +  I  f  (t)dx.  (63) 

Jn  Jon  Jn 

for  almost  all  t  G  (0 ,T)  where 

Q (t)  =  ^(Vsu h(t))  n  -  j^uh(t)^J 

for  a  given  mesh  face  F  ^Yh. 

By  interpreting  Nitsche’s  method  as  a  variational  wall  model  as  was  discussed  in 
Subsection  6.1,  we  see  that  the  above  balance  law  dictates  that  linear  momentum  en¬ 
ters  or  leaves  the  system  either  through  body  forces  or  surface  traction  forces.  Hence, 
our  semi-discrete  formulation  properly  mimics  the  continuous  problem.  We  can  derive 
similar  balance  laws  when  non-homogeneous  Dirichlet,  traction,  or  periodic  bound¬ 
ary  conditions  are  specified  instead  of  homogeneous  Dirichlet  boundary  conditions. 
Additionally,  we  can  derive  momentum  balance  laws  over  subdomains  by  introducing 
auxiliary  flux  spaces  on  subdomain  boundaries.  For  more  on  this  procedure,  see  [31]. 


7.2  Conservation  Properties  on  Cylindrical  Domains 

When  the  domain  fl  is  not  rectilinear,  our  semi-discrete  formulation  is  no  longer 
guaranteed  to  globally  conserve  linear  momentum.  This  is  because  the  Piola  trans¬ 
formation  does  not,  in  general,  map  constant  vector  fields  in  parametric  space  to 
constant  vector  fields  in  physical  space.  However,  the  special  structure  of  the  Pi¬ 
ola  transformation  allows  our  formulation  to  admit  more  natural  momentum  balance 
laws  for  general  parametric  domains.  In  this  subsection,  we  demonstrate  that  if  Ft  is 
a  cylindrical  domain  defined  via  a  polar  mapping,  our  formulation  globally  conserves 
axial  angular  momentum  and  axial  linear  momentum. 

To  proceed,  let  us  introduce  the  mapping 


F(6,6,6) 


{(rout  -  An)6  +  An)  sin(27r£i) 
((rout  ~  An)6  +  An)  cos(2vr£i) 
#6 


,V(6,6,6)e(0,i)3 


from  the  parametric  domain  H  =  (0,  l)3  to  a  physical  domain  H  set  between  two  con¬ 
centric  cylinders.  The  radius  of  the  inner  cylinder  is  taken  to  be  rm ,  the  radius  of  the 
outer  cylinder  is  taken  to  be  rout ,  and  the  heights  of  the  cylinders  are  taken  to  be  Ft. 
Periodicity  is  applied  in  the  direction.  We  have  the  following  relationship  between 
our  parametric  coordinate  system  and  the  cylindrical  coordinate  system  ( 9,r,z ): 


9  =  2tt6 

r  =  (r  out  -  An)  +  An 
*  =  ^  3. 
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Without  loss  of  generality,  let  us  assume  that  rm  =  1,  rout  =  2,  and  H  =  1.  The 
semi-discrete  velocity  functions  v/(,  G  14  have  the  representation 


Vfc(x(^))  :  = 


det  (£>F(£)) 

where  V/,  G  V/,,.  A  direct  computation  shows 


DF(i)  vfc(£) 


Vfc(x($))  := 


cos(d) 
—  sin(d) 
0 


0 


sin(fl) 

2nr 

cos(0)  n 

2n  r 

0  2V  -J 


Now,  since  the  space  14  contains  all  vector-valued  functions  which  are  linear  polyno¬ 
mials  in  £l,  £2,  and  £3  (since  k'  >  1),  we  have  thus  proven  the  vector- valued  functions 

COS(0)  fffl  0 

-sin  (9)  0 

0  0  A- 


s  = 


r 

r  cos(d) 

y 

0 

= 

— r  sin(d) 

= 

—x 

_  0  _ 

0 

0 

and 


z  = 


cos(d) 


sin  (8) 


2nr 

-MO) 

0  0  2W 


0 

0 


0 

"  0  " 

0 

= 

0 

27rr 

1 

are  members  of  V/,,.  Let  us  now  select  wh  =  s  and  qh  =  Phif)  in  (62).  Since 


Vs  = 


0  1  0 

-10  0 
0  0  0 


is  an  anti-symmetric  matrix,  we  have 


(2z/  (Vsu h(t))  :  (Vss)  -  (u h(t)  <g>  u h(t))  :  Vs)  dx 


Ida 


2v  ((Vss)n)  •  uh(t)ds  =  0 


for  almost  every  t  G  (0,T).  Therefore, 

4:  f  uh{t)  •  sdx  =  [  ( -Ptrace,h(t)n  +  2v  ( (VAu ,h(t))  n  -  y^uh(i) 
dt  Jn  Jon  V  V  hF 

+  [  f(f)  ■  sdx. 


■  sds 


Furthermore,  we  have  that 


w  ■  s  =  (w  x  r)  ,  Vw  G 
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where  (-)z  denotes  the  ^-component  of  a  vector  and 


r  = 


r  sin($) 
r  cos(d) 
0 


Hence,  we  have  proven  the  following  balance  law  for  axial  angular  momentum. 


Global  conservation  of  axial  angular  momentum  on  cylindrical  domains. 

For  cylindrical  domains,  the  semi-discrete  velocity  field  uh  satisfies 


i  f  K(i)xr)2rfx=  [  ((—ptrace,h(t)n  +  Q(t))  x  r),  ds 
(lt  Jn  Jon 

+  [  (fW  X  r)2rfx 
Jn 


(64) 


for  almost  all  t  e  (0,T)  where 


r  = 


r  sin(0) 
r  008(6*) 
0 


and 


=  2u  (Vsu/)(t))  n 


a 


pen 


hr 


U  hit) 


for  a  given  mesh  face  F  e  Th. 


The  above  balance  law  dictates  that  axial  angular  momentum  enters  or  leaves 
the  semi-discrete  system  either  through  applied  moments  or  torques.  This  balance 
law  mimics  the  corresponding  continuous  balance  law  for  axial  angular  momentum. 
We  would  like  to  mention  that  derivation  of  this  balance  law  was  contingent  upon 
employing  the  symmetrized  gradient  for  the  viscous  stress  term  instead  of  the  full 
gradient. 

Selecting  wh  =  z  and  =  Ph(t)  in  (62)  yields  the  following  balance  law  for  axial 
linear  momentum.  It  states  that  axial  momentum  enters  or  leaves  the  system  either 
through  axial  body  forces  or  axial  traction  forces. 


Global  conservation  of  axial  linear  momentum  on  cylindrical  domains.  For 

cylindrical  domains,  the  semi-discrete  velocity  field  uh  satisfies 

4:  /  (w(t))zdx=  [  ((-ptraceh(t)n  +  Q(t)))  ds 
dt  Jn  Jan 

+  [  (f0))2dx  (65) 

Jn 

for  almost  all  t  G  (0,T)  where 

Q (t)  =  2 v  ^(Vsu h(t))  n  -  h(t)^j 


24 


for  a  given  mesh  face  F  G  IV 

As  a  final  note,  we  would  like  to  state  that  if  we  had  utilized  a  NURBS  mapping 
to  represent  the  physical  domain  we  would  have  arrived  at  a  global  conservation 
statement  for  axial  angular  momentum  but  not  for  axial  linear  momentum,  ffence,  we 
believe  polar  mappings  hold  a  distinct  advantage  over  NURBS  mappings  when  solv¬ 
ing  problems  harboring  important  symmetry  and  conservation  properties.  This  being 
said,  we  would  like  to  mention  that  one  recovers  both  linear  and  angular  momentum 
balance  laws  when  employing  NURBS-based  isogeometric  analysis  in  conjunction  with 
a  residual-based  Variational  Multiscale  method  [3].  Such  a  discretization,  however, 
does  not  exactly  replicate  the  incompressibility  constraint  and  alternatively  attempts 
to  model  the  effects  of  fine-scale  solution  features  (including,  in  this  case,  compress¬ 
ibility)  on  the  resolved  components  of  the  flow. 


7.3  Vorticity,  Enstrophy,  and  Helicity 

We  continue  this  section  by  deriving  discrete  balance  laws  for  vorticity,  enstrophy, 
and  helicity  when  the  applied  forcing  term  is  conservative  (i.e.,  when  the  forcing  is 
derived  through  a  potential).  In  order  to  do  so,  we  must  properly  define  vorticity 
at  the  semi-discrete  level  (see,  for  example,  [43]).  We  restrict  our  discussion  to  the 
three-dimensional  setting  for  the  remainder  of  this  section.  Recall  that  the  vorticity 
held 

u;  =  curlu 

satisfies  the  partial  differential  equations 

duo 

—  +  V  •  (u  <g)  cj  —  u)  =  V  •  ( 2vWsoo )  (66) 

and 

diva;  =  0  (67) 

when  the  applied  forcing  term  is  conservative.  Furthermore,  since  u  =  0  on  <9h2, 

uj  ■  n  =  0.  (68) 


In  light  of  these  identities,  conservation  statements,  and  boundary  conditions,  we 
define  semi-discrete  vorticity  as  the  solution  of  the  following  problem:  find  a >h  G 
V}f  =  {v/j  G  Vf  :  divvfe  =  0}  such  that  o^(0)  =  a ;0,h  and 

(dtu:h(t),  vft)L2(n)  +  kh(u>h(t),Vh)  +  c(uh(t),ujh(t)]  vh )  -  c(ujh(t),  u h(t);  vh) 

-  £  L 2" 

F£  VhJF 

for  almost  every  t  G  (0,  T )  and  for  all  wh  G  Vo  h  where  u>0th  G  Vo  fc  is  a  suitably  defined 
initial  condition.  Note  that  we  have  strongly  enforced  normal  boundary  conditions  for 
the  semi-discrete  vorticity  held  and  weakly  enforced  tangential  boundary  conditions 


-  ((Vsvfe)  n)  •  curlu h(t)  +  curlu h(t)  ■  vh  )  ds  (69) 

flfP 
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using  Nitsche’s  method.  It  should  be  further  noted  that  the  semi-discrete  problem 
given  by  (69)  may  not  have  a  global-in-time  solution.  This  is  due  to  the  presence 
of  vortex  stretching.  However,  Picard’s  existence  theorem  can  be  used  to  show  the 
problem  has  a  unique  local-in-time  solution.  Hence,  let  us  assume  for  the  remainder 
of  this  discussion  that  the  end-time  T  has  been  sufficiently  chosen  such  that  the  semi¬ 
discrete  vorticity  equation  has  a  unique  solution.  Beyond  such  T,  the  semi-discrete 
vorticity  may,  in  theory,  experience  blow-up.  Such  a  hypothetical  blow-up  lies  at  the 
heart  of  the  Navier-Stokes  Millenium  Problem. 

We  now  present  a  discrete  balance  law  for  vorticity  which  holds  for  periodic  do¬ 
mains.  It  is  a  simple  consequence  of  choosing  =  e*  in  (69)  where  e*  is  a  unit 
vector.  Analogous  balance  laws  can  be  proven  for  rectilinear  domains  subject  to  no¬ 
penetration  and  no-slip  boundary  conditions. 


Global  balance  law  for  vorticity.  Suppose  fl  is  the  three-dimensional  torus.  For 
conservative  applied  forces,  the  semi-discrete  vorticity  field  u>h  satisfies 

d  f 

—  /  u>h(t)dx  =  0  (70) 

for  almost  all  t  G  (0,T). 


We  next  derive  a  discrete  balance  law  for  enstrophy.  To  begin,  we  insert  v/(  = 
u>h(t)  into  (69)  for  almost  every  t  G  (0,T): 

(0twft(t),u;ft(t))L2(n)  +  kh(ujh(t),ujh(t))  +  c(uh(t),u:h(t)-,ujh(t))  -  c(uh(t),uh(t);uh(t)) 


£ 

F&h' 


2  v 


((Vsu:h(t))  n)  •  curlu h(t)  +  ^^curluA(t)  •  u>h(t)  )  ds. 


A  straight-forward  calculation  invoking  the  divergence  theorem  and  taking  advantage 
of  the  fact  that  the  semi-discrete  velocity  field  is  divergence-free  yields 


(71) 


The  product  rule  gives 

d  f 

(dtu>h(t),u>h(t))L2(n)  =  -T7  7  h(t)dx. 

J  ^ 

where  7 hit)  =  -,\tjJh{t)\2  is  the  semi-discrete  enstrophy  density  of  the  fluid.  To  handle 
the  term  corresponding  to  vortex  stretching,  we  employ  integration  by  parts  and  the 
fact  that  diva;/!  =  0  to  write 


c(ujh(t),uh(t)-,ujh(t)) 
A  direct  calculation  then  gives 


u>h(t)TVu(t)ujh(t)dx. 


uh(t)TVu(t)ujh(t)  =  u:h(t)TB(uh(t))ujh(t) 
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where  B(u h(t))  =  Vsu/l(t)  is  the  semi-discrete  rate  of  strain  tensor.  Thus,  we  have 
arrived  at  the  following  balance  law. 


Global  balance  law  for  enstropy.  For  conservative  applied  forces,  the  semi¬ 
discrete  vorticity  field  u)h  satisfies 


f  lfh(t)dx  —  —  f  2v\Vsu>h(t)\2d-x  +  f  u:h(t)TI])(uh(t))ujh(t)dx 


+  /  2v  ((Wsu)h(t))  n)  ■  u>h(t)ds 
Jon 

a 


J2  [  2u  ((Vsw^(t))n 

Ferh  Jf  ' 


f^-uh{t)  )  ■  (ujh(t)  -  curluA(t))ds 

tip 


(72) 


for  almost  all  t  e  (0,T). 


The  first  and  second  lines  of  (72)  contain  enstrophy  production  terms  associated 
with  the  domain  interior  and  boundary  respectively  while  the  third  line  of  (72)  con¬ 
tains  supplemental  production  terms  associated  with  weak  enforcement  of  boundary 
conditions.  When  the  boundary  condition 

ioh  =  curlu/, 

is  met  exactly,  the  third  line  of  (72)  vanishes  and  we  recover  the  same  entropy  balance 
law  as  satisfied  by  the  exact  vorticity  field.  Note  that  since  the  divergence  of  U/(  is 
precisely  zero,  D(u^(f))  has  the  same  indefinite  structure  as  its  continuous  counter¬ 
part.  That  is,  it  has  three  real  eigenvalues  which  sum  to  zero.  Consequently,  we  have 
appropriately  captured  the  vortex  stretching  term  with  our  semi-discrete  formulation. 
The  same  cannot  be  said  for  semi-discretizations  which  satisfy  the  incompressibility 
constraint  only  approximately,  even  if  the  momentum  and  vorticity  equations  are 
written  in  skew-symmetric  form. 

We  finish  here  by  deriving  a  discrete  balance  law  for  helicity.  First  insert  (v^,  (p,)  = 
(i u>h(t),Ph(t ))  into  (40)  and  =  u hit)  into  (69)  for  almost  every  t  G  (0,T).  Adding 
the  two  resulting  expressions  and  taking  into  consideration  the  fact  that  f  is  conser¬ 
vative,  we  obtain 


(dtnh(t),  u?ft(f))L2(n)  +  kh(uh(t),u>h(t))  +  c(uh(t),  u h(t);  u }h(t)) 
+(dtuh(t),  uh(t))La(n)  +  kh(uh(t),  u h(t))  +  c(ufc(t),  u>h(t);  u h(t))  -  c(uh(t),  ufc(t);  u h(t)) 

=  J  2u  ((Vsu h(t))  n)  •  curlu h(t)  +  ^p^curlufc(t)  •  u h(t)^  ds. 


fg  rh 

By  the  product  rule,  we  have 


(dtuh(t),ujh(t))L2m  +  (0tu>A(t),uft(t))L2(n)  =  —  /  Qh(t)d 
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where  Qh(t)  =  uift(t)  -Uft(t)  is  the  semi-discrete  helical  density  of  the  flow.  Integration 
by  parts  and  the  fact  that  our  semi-discrete  velocity  field  is  divergence-free  give 

c(uh(t),uh(t)-,u:h(t))  =  -c(nh(t),u:h(t)-,uh(t)). 

Finally,  a  straight-forward  calculation  invoking  the  divergence  theorem  and  taking 
advantage  of  the  fact  that  the  semi-discrete  vorticity  field  is  divergence-free  yields 

c(ujh(t),uh(t);uh(t))  =  0. 

Collecting  our  equations,  we  acquire  the  following  balance  law  for  helicity. 


Global  balance  law  for  helicity.  For  conservative  applied  forces,  the  semi-discrete 
solution  satisfies 

^  [  Qh{t)dx  =  —  f  Av(\7su)h(t))  :  (Vsuh(f))dx 


idci 


Ida 


2z/((Vsu/l(f))n)  -uh(t)ds 
2 v  ((Vsu>h(t))  n)  •  uh(t)ds 


Y  I  2v  ({Vsuh(t))n- 

pery  Jf  V 

Y  f  2v  ({Vsuh(t))n- 

Pcr,  J F  k 


a 


pen 


hi 


uh(t)  ■  uh(t)ds 


a 


h(t)  )  ■  (u )h(t)  -  curlu h(t))ds 
hp 


(73) 


for  almost  all  t  e  (0,T). 


The  first  line  of  (73)  contains  helicity  production  terms  associated  with  the  domain 
interior,  the  second  and  third  lines  of  (73)  contain  production  terms  associated  with 
the  domain  boundary,  and  the  fourth  and  fifth  lines  of  (72)  contain  supplemental 
production  terms  associated  with  weak  enforcement  of  boundary  conditions.  When 
the  boundary  conditions 

uh  =  0 

and 

ujh  =  curl  u/, 

are  met  exactly,  the  fourth  and  fifth  lines  of  (73)  vanish  and  we  recover  the  same 
helicity  balance  law  as  satisfied  by  the  exact  vorticity  held. 

Note  that  in  the  limit  of  vanishing  viscosity  our  global  helicity  balance  reduces  to 

d  f 

It  =  o. 

Thus,  just  as  in  the  infinite-dimensional  setting,  helicity  is  an  inviscid  invariant  for 
our  semi-discrete  formulation.  We  believe  this  is  a  very  important  property  given  the 


pivotal  role  helicity  plays  in  flow  structure  development.  As  a  final  remark,  we  would 
like  to  mention  all  of  the  discrete  balance  laws  presented  here  generalize  to  other  sets 
of  boundary  conditions  including  non- homogeneous  Dirichlet  boundary  conditions, 
prescribed  traction  boundary  conditions,  and  periodic  boundary  conditions. 


8  Numerical  Experiments 

In  this  section,  we  numerically  test  our  semi-discretization  scheme  using  three  selected 
benchmark  problems:  two-dimensional  Taylor-Green  vortex  flow,  alternating  cylin¬ 
drical  Couette  flow,  and  three-dimensional  Taylor-Green  vortex  flow.  Throughout, 
we  choose  Nitsche’s  penalty  constant  as 

Cpen  =  5  (k‘  +  1) 

which  we  have  found  to  be  sufficiently  large  in  order  to  ensure  numerical  stability. 
Additionally,  we  employ  uniform  parametric  meshes  and  B-spline  spaces  of  maximal 
continuity. 


8.1  Two-Dimensional  Taylor-Green  Vortex  Flow 


As  a  first  numerical  experiment,  we  consider  two-dimensional  Taylor-Green  vortex 
flow.  Two-dimensional  Taylor-Green  vortex  flow  is  a  simple  periodic  (in  space)  vor¬ 
tical  flow  subject  to  the  initial  condition 


u0(a,2/) 


sin(x)  cos (y) 
—  cos(a;)  sin(y) 


The  exact  solution  for  this  flow  exponentially  decays  in  time  and  satisfies  the  rela¬ 
tionships 


u  {x,y,t) 
p(x,y,t) 


sin(x)  cos (y) 
—  cos  (a;)  sin(y) 


exp(— 2  ut), 


1 

4 


(cos(2a;)  +  cos(2y))  exp(— 4 ut). 


It  is  easily  seen  that  the  nonlinear  convection  term  is  exactly  balanced  by  the  pressure 
term  and  thus  does  not  interfere  with  the  evolution  of  the  velocity  flow  field.  Hence,  a 
question  of  practical  interest  is  whether  or  not  the  nonlinear  convection  term  interferes 
with  the  evolution  of  the  velocity  field  at  the  discrete  level. 

We  have  simulated  two-dimensional  Taylor- Green  vortex  flow  using  divergence- 
conforming  B-spline  discretizations  of  varying  mesh  size  and  polynomial  degree.  We 
restricted  our  computations  to  the  domain  =  (0, 7r)2  by  employing  symmetry  con¬ 
ditions  along  dfl.  A  linear  parametric  mapping  was  utilized  to  describe  the  physical 
domain.  The  Crank-Nicolson  method  [16]  was  employed  to  discretize  in  time,  and 
the  time-step  size  was  chosen  to  be 


At  =  min 


,  fc'+i 

h  2 


5 
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Table  1:  Two-dimensional  Taylor-Green  vortex  flow:  Convergence  rates  at  t  —  15  for 

f  kf  4-1  t  2  'l 

Re  =  100.  Time-step  size  chosen  as  At  =  min  <  h~^~ ,  ^  >. 


Polynomial  degree  k'  =  1 


h 

1/8 

1/16 

1/32 

1/64 

u  —  uh  HpQ) 

1.87e-l 

9.34e-2 

4.67e-2 

2.34e-2 

order 

- 

1.00 

1.00 

1.00 

u  —  uh  L2(0) 

1.02e-2 

2.51e-3 

6.24e-4 

1.56e-4 

order 

- 

2.02 

2.01 

2.00 

\\P-Ph\\v(n) 

1.09e-2 

2.57e-3 

6.32e-4 

1.59e-4 

order 

- 

2.08 

2.02 

1.99 

Polynomial  degree  k' 

=  2 

h 

1/8 

1/16 

1/32 

1/64 

u  —  Uh  H1(Q) 

9.64e-3 

2.38e-3 

5.92e-4 

1.48e-4 

order 

- 

2.02 

2.01 

2.00 

u  —  uh  L2(f2) 

5.96e-4 

7.24e-5 

8.98e-6 

1.12e-6 

order 

- 

3.04 

3.01 

3.00 

\\p-Ph\\v(n) 

1.39e-3 

1.56e-4 

1.90e-5 

2.36e-l 

order 

- 

3.16 

3.04 

3.01 

Polynomial  degree  k'  =  3 


h 

1/8  1/16  1/32 

U  —  Uh  H'ffi) 

order 

u  —  u/i  L2(Q) 

order 

\\p-Ph\\mn) 

order 

5.39e-4  6.88e-5  8.76e-6 

3.00  2.97 

3.42e-5  2.15e-6  1.36e-7 

3.99  3.98 

1.69e-4  9.44e-6  5.77e-7 
4.16  4.03 

30 


Figure  2:  Two-dimensional  Taylor-Green  vortex  flow:  L 2  error  of  velocity  field  versus 
time  at  Re  =  10,  20, 40,  80, 160, 320, 640, 1280,  oo  for  a  k'  —  1  discretization  with 
16  x  16  elements.  Time-step  size  chosen  as  At  =  min 


sufficiently  small  to  ensure  temporal  discretization  errors  are  of  the  same  order  as 
spatial  discretization  errors.  It  should  be  mentioned  the  use  of  the  Crank-Nicolson 
method  results  in  fully  discrete  schemes  which  automatically  inherit  the  conservation 
structure  of  the  corresponding  semi-discrete  schemes.  The  initial  condition  was  chosen 
using  L“-projection  into  the  discrete  space  of  divergence- free  velocity  fields.  Conver¬ 
gence  rates  obtained  at  time  t  =  15  for  a  flow  of  Reynold’s  number  Re  =  ^  =  100  are 
provided  in  Table  1.  Note  from  the  L2-  the  H1-norms  of  the  velocity  error  and  the 
L2-norm  of  the  pressure  error  optimally  converge  in  h.  To  analyze  the  behavior  of 
our  method  in  time,  we  have  plotted  the  L2-norm  of  the  velocity  error  versus  time  for 
a  chosen  spatial  discretization  and  for  a  wide  variety  of  Reynold’s  numbers  in  Figure 
2.  Note  from  the  figure  that  our  numerical  error  is  bounded  in  time.  Moreover,  the 
numerical  error  decays  roughly  at  the  same  rate  as  the  exact  solution.  Indeed,  we 
are  able  to  reproduce  a  time-indepedent  solution  when  Re  =  oo.  This  indicates  the 
nonlinear  convection  term  has  not  interfered  with  the  flow  evolution  of  our  discrete 
velocity  solution. 

To  contrast  our  methodology  with  standard  mixed  flow  discretizations,  we  re¬ 
peated  the  above  computations  for  a  conservative  Taylor- Hood  finite  element  approx¬ 
imation.  Again,  the  Crank-Nicolson  method  was  employed  to  discretize  in  time.  We 
found  that  the  results  obtained  using  this  flow  technology  were  unstable  in  general. 
To  illustrate  this,  we  have  plotted  in  Figure  3  the  L2-norm  of  the  velocity  error  ver¬ 
sus  time  for  a  Q2/Q1  Taylor- Hood  discretization  at  Re  =  00  on  a  mesh  with  8x8 
elements.  Note  the  exponential  blow-up  of  the  error  in  time.  This  blow-up  is  a  direct 
result  of  unphysical  energy  growth  stemming  from  the  nonlinear  convection  term.  In¬ 
deed,  we  have  been  unable  to  stably  compute  the  discrete  flow  solution  beyond  a  time 
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Figure  3:  Two-dimensional  Taylor-Green  vortex  flow:  Blow-up  of  the  L 2  velocity 
error  for  the  conservative  Q2/ Q\  Taylor- Hood  discretization  at  Re  —  oc  on  a  mesh 
with  8x8  elements.  Time-step  size  chosen  as  At  =  A, 

of  t  =  5.  These  results  are  a  testament  to  the  benefits  of  employing  a  conservative 
discretization  which  exactly  preserves  the  divergence- free  constraint. 

8.2  Alternating  Cylindrical  Couette  Flow 

As  a  second  numerical  experiment,  we  consider  the  flow  of  a  constant-property  New¬ 
tonian  fluid  lying  between  a  fixed  inner  cylinder  and  an  oscillating  outer  cylinder. 
This  flow  scenario  is  referred  to  as  alternating  cylindrical  Couette  flow.  The  problem 
setup  is  illustrated  in  Figure  4.  No  external  forcing  is  applied.  The  fluid  is  assumed 
to  be  at  rest  at  time  t  =  0.  Then,  the  outer  cylinder  begins  to  oscillate  with  angular 
velocity  equal  to  Uq  —  U sin  (cot),  inducing  the  fluid  to  slip  along  with  the  outer  cylin¬ 
der.  As  time  evolves,  the  flow  held  throughout  the  region  between  the  two  cylinders 
approaches  a  periodic  (in  time)  steady  state.  The  how  velocity  associated  with  this 
steady  state  can  be  explicitly  derived  (see,  for  example,  [49])  and  is  equivalent  to 

ug(r,t)  sin(6l) 

U  ug(r,t)  cos(6) 


with 


ug(r,  t )  =  U Irnag 


loMKoftrjn)  -  Io('yrin)K0('yr) 

logout) Kobrin)  -  I0 (qr in) K0 (jr ^ ) 


exp  {icot} 


where  (r,  6)  are  polar  coordinates  with  respect  to  the  center  of  the  cylinders,  7  = 
\Jicov,  and  Iq  and  K0  are  modified  Bessel  functions  of  the  hrst  and  second  kind, 
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Figure  4:  Alternating  cylindrical  Couette  flow:  Problem  setup. 

respectively.  Unfortunately,  no  closed-form  solution  exists  for  the  pressure  field.  The 
Reynold’s  number  for  this  flow  is  taken  to  be 


Re  = 


In  what  follows,  we  assume  rm  =  1,  Tout  =  2,  and  U  —  1. 

In  Figure  5,  we  have  plotted  at  time  instances  t  =  40.0,  40.2,  40.4,  40.6,  40.8, 
41.0  the  exact  angular  velocity  field  associated  with  a  Re  =  200  flow  subject  to 
an  oscillation  frequency  of  u  =  5.  At  these  particular  time  instances,  the  flow  has 
already  reached  the  steady  periodic  state.  Note  from  the  figure  that  there  is  a  small 
boundary  layer  attached  to  the  outer  cylinder.  Further  note  that  there  is  substantial 
flow  reversal  in  a  region  away  from  the  outer  cylinder. 

We  believe  alternating  cylindrical  Couette  flow  is  an  interesting  and  challenging 
numerical  test  problem  for  a  number  of  reasons.  First  of  all,  the  problem  exhibits 
important  symmetries  that  ideally  should  be  preserved  in  a  numerical  simulation.  As 
a  consequence  of  these  symmetries,  the  nonlinear  advection  term  is  exactly  balanced 
by  the  pressure  term.  Second,  the  problem  is  characterized  by  strong  shifts  in  angular 
momentum  in  time.  Consequently,  a  methodology  which  admits  angular  momentum 


33 


Figure  5:  Alternating  Couette  flow:  Plot  of  the  angular  velocity  field  at  Re  =  200  for 
t  =  40.0, 40.2, 40.4, 40.6, 40.8, 41.0. 
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Table  2:  Alternating  cylindrical  Couette  flow:  Convergence  rates  at  t  —  40  for  Re 

fc/  +  l 

200  and  uj  —  5.  Time-step  size  chosen  as  At  =  h^r . 


Polynomial  degree  k!  =  1 


h/h  o 

1/8 

1/16 

1/32 

1/64 

u  —  1  H'(Q) 

2.99e0 

1.79e0 

1.05e0 

5.51e-l 

order 

- 

0.74 

0.77 

0.93 

11  —  uh  1  L2(0) 

9.35e-2 

2.95e-2 

8.86e-3 

2.40e-3 

order 

- 

1.66 

1.74 

1.88 

||  Ur  i^r)h  ||  L2(f2) 

0 

0 

0 

0 

Polynomial  degree  k' 

=  2 

h/h0 

1/8 

1/16 

1/32 

1/64 

u  —  HflQ) 

2.24e0 

8.72e-l 

2.14e-l 

5.15e-2 

order 

- 

1.36 

2.03 

2.06 

u  —  uh  |  L2(0) 

3.23e-2 

7.00e-3 

9.56e-4 

1.26e-4 

order 

- 

2.21 

2.87 

2.92 

\\ur  —  (ur)h\  L2(Q) 

0 

0 

0 

0 

Polynomial  degree  h!  =  3 


h/h0 

1/8 

1/16 

1/32 

U  —  Uh  Hi(f2) 

1.51e0 

1.84e-l 

2.32e-l 

order 

- 

3.04 

2.99 

u  —  u/i  L2(Q) 

2.13e-2 

1.27e-3 

7.94e-5 

order 

- 

4.07 

4.00 

\\ur  —  (ur)h\  L2(H) 

0 

0 

0 
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Figure  6:  Alternating  cylindrical  Couette  flow:  L2  error  of  velocity  field  versus  time 
at  Re  =  200  for  a  k!  =  1  discretization  with  128  x  32  elements  in  the  6 ,  r  directions. 
The  oscillation  frequency  for  this  simulation  was  chosen  as  cu  =  5,  and  the  time-step 
size  was  chosen  as  At  —  h. 

balance  is  preferred.  Finally,  the  problem  is  characterized  by  the  presence  of  boundary 
layers  and  flow  reversal.  Many  flow  technologies  exhibit  spurious  oscillations  when 
applied  to  problems  harboring  such  features. 

We  have  numerically  simulated  alternating  cylindrical  Couette  flow  using  divergence- 
conforming  B-spline  discretizations  of  varying  mesh  size  and  polynomial  degree.  We 
utilized  a  polar  mapping  to  represent  the  annular  domain,  and  Nitsche’s  method  was 
invoked  to  enforce  the  slip  condition  along  the  cylinder  surfaces.  The  resulting  semi¬ 
discretizations  satisfy  an  angular  momentum  balance  law  as  discussed  in  Subsection 
7.2.  As  in  the  last  verification  test,  the  Crank-Nicolson  method  was  employed  to 
discretize  in  time,  and  the  time-step  size  was  chosen  to  be 


sufficiently  small  to  ensure  temporal  discretization  errors  are  of  the  same  order  as 
spatial  errors.  Convergence  rates  obtained  at  t  =  40  for  a  Re  —  200  flow  subject  to 
an  oscillation  frequency  of  u  =  5  are  provided  in  Table  2.  Note  from  the  table  that  the 
LJ-  and  H1-norms  of  the  velocity  error  approach  optimal  convergence  rates  in  h  and 
that  we  have  obtained  axisymmetric  velocity  fields  with  zero  radial  component.  To 
analyze  the  behavior  of  our  method  in  time,  we  have  plotted  in  Figure  6  the  L2-norm 
of  the  velocity  error  versus  time  for  a  chosen  spatial  discretization.  Note  from  the 
figure  that  our  numerical  error  is  bounded  and  periodic  in  time.  This  indicates  that 
our  numerical  solution,  like  the  exact  solution,  has  reached  a  periodic  steady-state. 
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8.3  Three-Dimensional  Taylor-Green  Vortex  Flow 

As  a  third  and  final  numerical  experiment,  we  consider  three-dimensional  Taylor- 
Green  vortex  flow.  Three-dimensional  Taylor-Green  vortex  flow  is  one  of  the  simplest 
systems  in  which  one  can  study  enstrophy  production  and  the  turbulence  resulting 
from  vortex  stretching.  The  initial  conditions  for  this  flow  are 


uo  (x,y,z) 


sin(x)  cos (y)  cos(T) 
—  cos(T)  sin (y)  cos(z) 


0 


The  flow  is  periodic  in  all  three  spatial  directions  in  the  domain  hi  =  (0,  2i r)3  and 
exhibits  a  64-fold  symmetry  which  can  be  exploited  in  numerical  simulation  [9].  The 
Reynolds  number  for  this  flow  is  commonly  taken  to  be 


Re  = 


1 

v 


In  Figure  7,  we  have  reproduced  time  history  plots  of  the  dissipation  rate 


e  =  — yr  /  2i/|Vsu|2dx  =  —  /  z/|Vu|2dx=  —  /  v\ curlu|2dx 


m 


i«i 


M 


that  were  obtained  by  Brachet  et  al.  in  [9]  via  Fourier-based  Direct  Numerical  Sim¬ 
ulation  (DNS)  with  2563  resolved  modes.  Note  that  the  flow  exhibits  significant 
enstrophy  production  throughout  the  initial  stages  of  flow  evolution  regardless  of 
Reynold’s  number.  At  Re  =  100,  the  time  corresponding  to  the  maximum  dissipa¬ 
tion  rate  is  approximately  t  ~  4.75.  As  the  Reynolds  number  is  increased,  the  time 
corresponding  to  the  maximum  dissipation  rate  gradually  increases  until  it  settles 
around  a  value  of  t  ~  9. 

To  simulate  three-dimensional  Taylor-Green  vortex  flow,  we  have  utilized  divergence- 
free  B-spline  discretizations  of  varying  mesh  size  and  polynomial  degrees  k'  =  1,  2,  3. 
We  have  exploited  symmetry  conditions  in  order  to  solve  the  unsteady  Navier-Stokes 
equations  on  the  restricted  domain  (0,  n)2  and  hence  reduce  the  dimensionality  of  our 
discrete  system  by  a  factor  of  8.  A  linear  parametric  mapping  was  utilized  to  describe 
the  physical  domain.  The  Crank-Nicolson  method  was  employed  to  discretize  viscous 
terms  while  the  Adams-Bashforth  multi-step  method  [2]  was  employed  to  discretize 
the  nonlinear  convective  terms.  This  time  discretization  procedure  allowed  us  to  avoid 
the  solution  of  a  nonlinear  algebraic  problem  at  each  time-step  at  the  cost  of  losing 
some  of  the  conservation  structure  of  the  semi-discrete  method.  A  time-step  size  of 
At  =  0.05/r  was  employed  in  all  of  our  simulations.  The  initial  condition  was  selected 
using  L2-projection  into  the  discrete  space  of  divergence-free  velocity  fields. 

In  Figure  8,  we  have  depicted  an  enstrophy  isosurface  associated  with  the  initial 
condition,  and  in  Figure  9,  we  have  depicted  an  enstrophy  isosurface  that  was  obtained 
at  time  t  =  6  via  a  third-order  B-spline  simulation  of  Re  =  200  flow  on  a  spatial  mesh 
comprised  of  32  x  32  elements.  This  time  roughly  corresponds  to  the  moment  of 
maximum  dissipation  rate.  Note  from  the  figures  that  while  the  initial  solution  is 
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Figure  7:  3-D  Taylor-Green  vortex  flow:  Time  history  plots  of  the  dissipation  rate  for 
various  Reynold’s  numbers.  Image  reproduced  from  Brachet  et  al.  [9]  with  permission 
from  Cambridge  University  Press. 
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Figure  8:  3-D  Taylor-Green  vortex  flow:  Visualization  of  enstrophy  isosurface  colored 
by  vertical  vorticity  at  t  —  0  for  Re  =  200.  (a)  3-D  View,  (b)  Overhead  view. 

Visualization  is  restricted  to  the  domain  (0, 7 r)3. 
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Figure  9:  3-D  Taylor-Green  vortex  flow:  Visualization  of  enstrophy  isosurface  colored 
by  vertical  vorticity  at  t  —  6  for  Re  =  200.  (a)  3-D  View,  (b)  Overhead  view. 

Visualization  is  restricted  to  the  domain  (0, 7 r)3. 
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(b) 


Figure  10:  3-D  Taylor-Green  vortex  flow:  Convergence  of  dissipation  rate  time  histo¬ 
ries  for  Re  =  200.  (a)  Convergence  of  k'  =  1  discretizations  under  mesh  refinement, 
(b)  Convergence  of  h  =  1/32  discretizations  under  degree  elevation. 
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comprised  of  a  single  vortex  on  the  restricted  domain  (0, 7r)3,  vortex  stretching  has 
separated  the  initial  vortex  into  many  vortical  structures  by  time  t  —  6.  Further  note 
the  vast  amount  of  symmetry  exhibited  by  the  vortical  structures.  We  found  that 
this  symmetry  was  preserved  in  all  of  our  numerical  experiments.  In  Figure  10(a),  we 
have  depicted  the  Re  =  200  dissipation  rate  time  histories  associated  with  a  sequence 
of  refined  k'  =  1  discretizations.  The  dissipation  rate  time  history  on  the  finest  mesh 
is  virtually  indistinguishable  from  the  corresponding  DNS  time  history  depicted  in 
Figure  7.  The  other  dissipation  rate  time  histories  quickly  converge  in  h.  It  should  be 
noted  that  we  have  been  able  to  stably  compute  arbitrary  Reynold’s  number  flow  on 
the  coarse  mesh  ( h  =  1/16),  though  the  results  were  wildly  inaccurate  at  long  times 
due  to  a  fine-scale  pile-up  of  energy  resulting  from  a  lack  of  resolution.  In  Figure 
10(b),  we  have  depicted  the  Re  =  200  dissipation  rate  time  histories  associated  with 
h  =  1/32  discretizations  of  varying  polynomial  degree.  Note  that  the  dissipation  rate 
time  histories  quickly  converge  in  k' .  Furthermore,  the  k'  =  3  dissipation  rate  time 
history  nearly  matches  the  corresponding  DNS  time  history  illustrated  in  Figure  7, 
though  the  k'  —  1,  h  —  1/64  results  are  in  slightly  better  agreement. 


9  Conclusions 

In  this  paper  we  have  described  the  use  of  divergence-conforming  B-spline  discretiza¬ 
tions  for  unsteady  Navier-Stokes  flows.  These  functions  enable  smooth,  pointwise 
divergence-free  solutions  to  be  computed  on  geometrically  mapped  meshes.  Conse¬ 
quently,  conservation  of  mass  is  satisfied  exactly,  both  globally  and  locally.  Our  semi¬ 
discrete  variational  equations  are  written  in  conservation  form  and  thus  preserve  other 
important  conservation  and  balance  laws.  This  is  in  fact  a  consequence  of  the  point- 
wise  exactness  of  mass  conservation.  We  note  that  heretofore  no  numerical  method 
for  the  Navier-Stokes  equations  has  been  able  to  achieve  similar  attributes.  Our  fo¬ 
cus  in  the  work  herein  has  been  on  Galerkin  discretizations  with  weakly-enforced 
no-slip  boundary  conditions.  To  be  more  precise,  we  satisfy  no-penetration  boundary 
conditions  strongly  and  the  tangential  boundary  condition  weakly.  As  we  have  docu¬ 
mented  in  previous  works,  this  treatment  of  the  no-slip  boundary  condition  is  always 
at  least  as  accurate  as  strong  enforcement  and,  in  the  case  of  under-resolved  bound¬ 
ary  layer  phenomena,  often  procedures  remarkably  accurate  results.  Furthermore, 
this  approach  provides  a  natural  pathway  to  the  computation  of  Euler  flows  as  the 
kinematic  viscosity  takes  on  the  limiting  value  of  zero.  In  the  computational  setting, 
no  changes  to  the  degree-of-freedom  structures  need  to  be  made  when  calculating 
Navier-Stokes  and  Euler  flows,  which  proves  very  convenient. 

The  semi-discrete  equations  are  shown  to  conserve  linear  and  angular  momen¬ 
tum,  and  discrete  balance  laws  for  vorticity,  enstrophy,  and  hclicity  are  also  derived. 
The  geometrical  structure  of  solutions  of  the  Navier-Stokes  and  Euler  equations  are 
intimately  linked  to  the  preservation  of  the  these  conservation  and  balance  laws. 
To  illustrate  the  behavior  of  the  methods  on  solutions  in  which  geometrical  struc¬ 
ture  is  sensitive  to  the  evolution  of  such  quantities,  we  performed  numerical  calcula¬ 
tions  of  two-dimensional  Taylor-Green  vortical  flow,  alternating  cylindrical  Couette 
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flow,  and  three-dimensional  Taylor-Green  vortical  flow.  In  each  case  the  structure¬ 
preserving  behavior  and  quantitative  accuracy  are  clearly  evident,  even  on  relatively 
coarse  meshes.  In  addition  to  the  computational  verifications  of  the  methodology, 
there  are  important  theoretical  implications:  this  is  potentially  the  first  instance  of 
a  structure-preserving  methodology  that  might  be  used  for  DNS  calculations  on  geo¬ 
metrically  interesting  domains.  As  such,  it  provides  a  numerical  conduit  for  exploring 
fundamental  flow  physics  on  geometries  never  before  explored. 

This  is  a  first  but  important  step  in  a  new  area  of  computational  fluid  dynam¬ 
ics.  To  go  further,  much  basic  research  needs  to  be  pursued.  To  develop  LES-level 
methods,  stabilized  and  variational  multiscale  generalizations  need  to  be  developed, 
preserving  the  conservation/balance  law  structure,  but  extending  it  to  arbitrarily  high 
Reynolds  numbers  on  mesh  sizes  that  are  approachable  with  contemporary  computer 
hardware.  In  addition,  the  methods  need  to  be  made  applicable  to  unstructured 
meshes  with  T-junctions  (i.e.,  “T-meshes”)  and  “extraordinary  points”.  We  hope  to 
play  a  significant  role  in  these  important  pursuits. 
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